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We  present  a  unifying  theoretical  framework  that  describes  recently  observed  many-body  effects  during  the 
interrogation  of  an  optical  lattice  clock  operated  with  thousands  of  fermionic  alkaline  earth  atoms.  The  frame¬ 
work  is  based  on  a  many-body  master  equation  that  accounts  for  the  interplay  between  elastic  and  inelastic 
p-wave  and  s-wave  interactions,  finite  temperature  effects  and  excitation  inhomogeneity  during  the  quantum 
dynamics  of  the  interrogated  atoms.  Solutions  of  the  master  equation  in  different  parameter  regimes  are  pre¬ 
sented  and  compared.  It  is  shown  that  a  general  solution  can  be  obtained  by  using  the  so  called  Truncated 
Wigner  Approximation  which  is  applied  in  our  case  in  the  context  of  an  open  quantum  system.  We  use  the 
developed  framework  to  model  the  density  shift  and  decay  of  the  fringes  observed  during  Ramsey  spectroscopy 
in  the  JILA  8'  Sr  and  NIST  171  Yb  optical  lattice  clocks.  The  developed  framework  opens  a  suitable  path  for 
dealing  with  a  variety  of  strongly-correlated  and  driven  open-quantum  spin  systems. 
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I.  INTRODUCTION 


One  of  the  ultimate  goals  of  modem  physics  is  to  understand  and  fully  control  quantum  mechanical  systems  and  to  exploit 
them  both,  at  the  level  of  basic  research  and  for  numerous  technological  applications  including  navigation,  communications, 
network  management,  etc.  To  accomplish  these  objectives,  we  aim  at  developing  the  most  advanced  and  novel  measurement 
techniques  capable  of  probing  quantum  matter  at  the  fundamental  level. 

Some  years  ago,  the  second  -  the  international  unit  of  time  -  was  defined  by  the  Earth’s  rotation.  However,  with  the  discovery 
of  quantum  mechanics  and  the  quantized  nature  of  the  atomic  energy  levels,  it  became  clear  that  atomic  clocks  could  be  more 
accurate  and  more  precise  than  any  mechanical  or  celestial  reference  previously  known  to  man.  Thus,  in  1967  the  second  was 
redefined  as  the  duration  of  9,192,631,770  periods  of  the  radiation  corresponding  to  the  transition  between  the  two  hyperfine 
energy  levels  of  a  caesium  atom.  Since  then,  the  accuracy  of  atomic  clocks  has  improved  dramatically,  by  a  factor  of  10  or  so 
every  decade.  The  characterization  of  the  unit  of  time  plays  a  central  role  within  the  International  System  of  Units  (SI)  because 
of  its  unprecedented  high  accuracy  and  because  it  is  also  used  in  the  definitions  of  other  units  such  as  meter,  volt  and  ampere. 

Thanks  to  the  development  of  laser  trapping  and  cooling  techniques  [1,2],  the  best  caesium  standards  have  reached  an  accuracy 
of  one  part  in  1016.  However,  caesium  clocks  are  limited  by  the  fact  that  they  are  based  on  atomic  transitions  in  the  microwave 
domain.  Because  the  quality  factor  of  the  clock  is  proportional  to  the  frequency,  optical  clocks  with  frequencies  that  can  be  106 
times  higher  than  microwaves,  offer  an  impressive  potential  gain  over  their  microwave  counterparts.  Optical  frequencies  on  the 
other  hand  are  very  difficult  to  measure,  as  the  oscillations  are  orders  of  magnitude  faster  than  what  electronics  can  measure.  The 
implementation  of  frequency  comb  technology  [3]  has  provided  a  coherent  link  between  the  optical  and  microwave  regions  of 
the  electromagnetic  spectrum,  greatly  simplifying  optical  frequency  measurements  of  high  accuracy.  After  the  development  of 
frequency  combs,  the  interest  in  optical  clocks  has  grown  rapidly.  Now,  optical  clocks  based  on  single  trapped  ions  and  neutral 
atoms  are  the  new  generation  of  frequency  standards  with  a  sensitivity  and  accuracy  as  high  as  one  part  in  1018  [4—6]. 

Optical  clocks  operated  with  fermionic  neutral  alkaline  earth  atoms  (AEA),  such  as  8 '  Sr  or  171  Yb,  have  matured  considerably. 
Those  employ  an  optical  lattice  to  tightly  confine  the  atoms  so  that  Doppler  and  photon-recoil  related  effects  on  the  transition 
frequency  are  eliminated.  State-of-the-art  neutral-atom-optical  clocks  have  surpassed  the  accuracy  of  the  Cs  standard  [7]  and 
just  recently,  thanks  to  advances  in  modern  precision  laser  spectroscopy,  are  reaching  and  even  surpassing  the  accuracy  of  single 
ion  standards  [6].  The  most  stable  of  these  clocks  now  operate  near  the  quantum  noise  limit  [8,  9].  The  stability  arises  from  the 
intrinsic  atomic  physics  of  two-valence-electron  atoms  that  possess  extremely  long  lived  singlet  and  triplet  states  (clock  states), 
with  intercombination  lines  nine  orders  of  magnitude  narrower  than  a  typical  dipole-allowed  electronic  transition. 

The  potential  advantage  of  neutral-atom  clocks  over  single  trapped  ion  clocks  is  that,  in  the  former,  a  large  number  of  atoms 
is  simultaneously  interrogated.  This  could  lead  to  a  large  signal-to-noise  improvement;  however,  high  atom  numbers  combined 
with  tight  confinement  also  lead  to  high  atomic  densities  and  the  potential  for  non-zero  collisional  frequency  shifts  via  contact 
atom-  atom  interactions.  With  atom-light  coherence  times  reaching  several  seconds,  even  very  weak  interactions  (e.g.,  fractional 
energy  level  shifts  of  order  >  1  x  10-16)  can  dominate  the  dynamics  of  these  systems. 

To  suppress  these  interactions,  the  use  of  ultracold,  spin-polarized  fermions  was  proposed.  The  idea  was  to  exploit  the 
Fermi  suppression  of  s-wave  contact  interactions  while  freezing  out  p-wave  and  higher  wave  collisions  at  ultracold  atomic 
temperatures.  Indeed,  at  precision  level  of  10~15,  the  JILA  Sr  clock  did  not  exhibit  a  density-dependent  frequency  shift  [10], 
however  as  the  measurement  precision  progressed,  density-dependent  frequency  shifts  were  measured  in  spin  polarized  fermions, 
at  the  JILA  Sr  clock  [1 1-13]  and  at  the  NIST  Yb  clock  [14], 

When  those  density-dependent  frequency  shifts  were  first  observed,  they  were  attributed  to  s-wave  collisions  allowed  by  inho¬ 
mogeneous  excitation  [11,  15-18],  under  the  assumption  that  p-wave  interactions  were  suppressed  at  the  operating  temperatures 
(T  ~  pK)  [11],  The  basic  initial  understanding,  obtained  from  a  mean-field  treatment,  was  that  excitation  inhomogeneities 
induced  by  the  optical  probing  laser  made  the  initially  indistinguishable  fermionic  atoms  distinguishable,  and  thus  allowed  them 
to  interact  via  s-wave  collisions. 

However,  studies  of  the  cold  collision  shift  in  the  NIST  Yb  optical  lattice  clock  using  Ramsey  spectroscopy  revealed  that 
p- wave  interactions  were  the  dominant  elastic  interactions  in  that  system  [19].  Furthermore,  evidence  of  inelastic  p-wave  inter¬ 
actions  was  reported  in  both  Yb  and  Sr  atomic  clocks  [20,  21].  Although  the  importance  of  many-body  interactions  in  optical 
clocks  has  been  recognized  theoretically  [16-18],  only  recent  measurements  have  revealed  their  many-body  nature  [22,  23].  In 
those  measurements,  the  role  of  s-wave  collisions  was  further  suppressed  by  operating  the  Sr  clock  with  highly  homogeneous 
atom-laser  coupling.  This  results  in  dominant  p- wave  interactions  with  a  collective  character,  as  we  will  explain  below. 

At  this  point  it  is  important  to  emphasize  that  recent  advances  in  modern  precision  laser  spectroscopy,  with  record  levels 
of  stability  and  residual  laser  drift  less  than  mHz/s  [8,  23,  24]  are  the  crucial  developments  that  are  allowing  us  to  deal  with 
clocks  operated  at  a  very  different  conditions  than  those  ones  dealt  with  just  few  years  ago.  The  level  of  precision  spectroscopy 
achievable  in  current  atomic  clocks  is  now  providing  the  required  spectral  resolution  to  systematically  resolve  and  study  the 
complex  excitation  spectrum  of  an  interacting  many-body  system.  This  was  certainly  not  the  case  in  prior  clock  experiments 
where  interaction  effects  were  subdominant  and  where  a  mean-field  treatment  was  more  than  enough  to  describe  the  clock 
behavior.  For  example  in  2006,  a  2  Hz  spectral  resolution  has  achieved  for  the  Sr  atomic  transition  and  no  interaction  effects 
were  observable  at  the  time  [25]. 


3 


FIG.  1.  (Color  online)  General  design  of  an  optical  lattice  clock 


In  this  paper  we  present  a  unifying  theoretical  framework  that  goes  beyond  a  simple  mean-field  treatment  and  that  is  capable 
of  describing  the  full  many-body  dynamics  of  nuclear  spin-polarized  alkaline  earth  atoms  during  clock  interrogation.  The  two 
clock  states  are  treated  as  an  effective  spin  degree  of  freedom.  Both,  elastic  and  inelastic  two-body  collisions  and  single-particle 
losses  are  present  during  the  dynamics,  and  thus  a  pure  Hamiltonian  formulation  is  not  sufficient.  Instead,  we  develop  a  master 
equation  formulation  which  is  capable  of  treating  the  quantum  evolution  of  an  open  spin  system.  We  provide  analytic/exact 
solutions  of  the  master  equation  dynamics  in  parameter  regimes  where  exact  treatments  are  possible.  For  the  more  generic 
situations  we  solve  the  dynamics  relying  on  the  so  called  Truncated  Wigner  Approximation  (TWA)  [26,  27],  In  contrast  to 
previous  theoretical  treatments  of  the  clock  dynamics,  which  were  limited  to  treating  two-particles  or  many-particles  but  at  the 
mean-held  level  or  under  the  all-to-all  approximation  [13,  16-18],  the  TWA  method  allows  us  to  include  both  elastic  and  inelastic 
collisions  beyond  the  mean-held  level,  hnite  temperature  effects  and  inhomogeneities  generated  by  either  the  laser  during  the 
pulse  interrogation  or  by  many-body  interactions.  Those  are  shown  to  be  crucial  for  properly  modeling  observed  many-body 
dynamics,  especially  at  T  >  10//.K,  at  which  excitation  inhomogeneities  can  not  be  neglected.  To  our  knowledge  this  is  the  hrst 
time  that  the  TWA  is  applied  to  describe  an  open  quantum  system  in  the  presence  of  inelastic  losses. 

Although  this  paper  focuses  on  optical  lattice  clocks,  the  developed  theoretical  framework  is  generic  for  driven  open-quantum 
systems  and  should  be  a  useful  platform  for  dealing  with  a  variety  of  current  experimentally  relevant  systems  including  trapped 
ions  [28,  29],  polar  molecules  [30-35],  nitrogen  vacancy  centers  [36],  and  atoms  in  optical  cavities  [37,  38]  among  others. 

The  remainder  of  the  paper  is  organized  as  follows.  In  section  II,  we  introduce  the  reader  to  the  basic  operation  of  an  atomic 
clock  and  derive  the  many-body  Hamiltonian  that  describes  the  dynamics  of  nuclear  spin-polarized  fermionic  atoms  during 
clock  interrogation.  We  then  proceed  to  derive  a  simplified  effective  spin  model  which  assumes  frozen  motional-excitations 
during  the  dynamics.  In  section  III,  we  solve  for  the  dynamics  under  the  assumption  of  collective  spin  interactions  (all-to-all 
interactions)  in  a  closed  system.  In  Section  IV,  we  show  how  to  treat  the  observed  two-body  losses  and  introduce  a  master 
equation,  which  we  solve  under  the  collective-interactions  approximation.  We  also  show  how  to  use  the  TWA  to  deal  with 
the  open  quantum  system  dynamics.  In  section  V,  we  relax  the  frozen-motional-degrees-of-freedom  approximation  and  derive 
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an  improved  spin  Hamiltonian  with  cubic  spin-spin  interactions  which  account  for  the  virtual  occupation  of  excited  motional 
modes.  In  section  VI,  we  go  beyond  the  all-to-all  interactions  approximation  and  present  a  more  general  prescription  that  can 
address  both  non-collective  interactions  and  losses  and  single-particle  inhomogeneities.  In  section  VII,  we  apply  the  developed 
theoretical  framework  to  model  the  dynamics  during  Ramsey  spectroscopy  observed  in  the  JILA8 '  Sr  and  the  NIST  171  Yb  optical 
lattice  clocks,  and  finally,  in  section  VIII,  we  present  the  conclusions.  In  Appendixes  1-5,  we  present  some  details  omitted  in  the 
main  text. 


II.  MANY-BODY  PHYSICS  DURING  CLOCK  OPERATION 
A.  A  simple  overview  of  an  optical  lattice  clock 

The  general  design  of  an  optical  lattice  clock  is  shown  in  Fig.  1.  It  consists  of  two  components,  a  laboratory  radiation  source 
and  an  atomic  system  with  a  natural  reference  frequency  determined  by  quantum  mechanics  to  which  the  laboratory  radiation 
source  can  be  compared.  Here,  the  laboratory  radiation  source  is  an  ultra-stable  continuous-wave  laser.  It  acts  as  the  local 
oscillator  (or  pendulum)  for  the  clock  and  is  used  to  probe  an  electromagnetic  resonance  in  an  atom.  The  atomic  signal  can  then 
be  used  to  determine  the  difference  between  the  laser  frequency  and  that  of  the  reference  atom,  allowing  laser  frequency  to  be 
monitored  and  stabilized  to  the  preferred  value.  A  frequency  comb  [3]  provides  the  gears  of  the  clock,  allowing  measurements 
of  the  laser  frequency  relative  to  other  high  accuracy  clocks  in  either  the  optical  or  microwave  domains. 

The  two  main  quantities  that  characterize  the  performance  of  a  clock  are  the  accuracy  and  the  precision.  The  accuracy 
is  determined  by  how  well  the  measured  frequency  matches  that  of  the  atom’s  natural  frequency.  In  general,  the  accuracy  will 
depend  on  the  atomic  species  used  and  how  well  it  can  be  isolated  from  environmental  effects  during  spectroscopy.  The  precision 
of  the  clock  is  more  commonly  referred  to  as  the  stability,  which  represents  the  repeatability  of  the  measured  clock  frequency 
over  a  given  averaging  time  r.  For  quantum  projection  noise  limited  measurements,  it  is  typically  expressed  as  [39] 


Here,  =  Q  is  the  line  quality  factor  of  the  clock  transition  for  a  linewidth  Au  and  \/N  is  associated  with  the  signal-to 
noise-ratio  achieved  for  interrogating  N  atoms  in  the  measurement  cycle  time  tc. 

Fermionic  AEA  such  as  8 'Sr  and  171  Yb  have  unique  properties  that  make  them  ideal  candidates  for  the  realization  of  atomic 
clocks  [40],  The  clock  states  are  the  ground  singlet  state,  1  So,  and  a  long  lived  triplet  state  3P0,  with  intercombination  lines 
both  electric  and  magnetic  dipole  forbidden  and  as  narrow  as  a  few  mHz,  see  Fig.  2.  In  the  ground  state  (15o),  and  to  leading 
order  in  the  excited  state  (3Pq),  the  electronic  degrees  of  freedom  have  neither  spin  nor  orbital  angular  momentum  [25],  This 
means  that  the  atoms  in  the  clock  states  only  interact  with  external  magnetic  fields  through  the  nuclear  spin  degrees  of  freedom 
which  have  a  (/-factor  1000  times  smaller  than  the  electronic  orbital  one.  Therefore,  AEA  are  much  less  sensitive  to  magnetic 
field  fluctuations  and/or  to  intensity  and  phase  noise  on  the  optical  fields  than  conventional  alkali  atoms. 

The  clock  transition  is  only  allowed  (i.e.  laser  light  weakly  couples  3So  to  3  Po)  because  in  the  excited  state  the  hyperfine 
interaction  leads  to  a  small  admixture  of  the  higher-lying  P  states  [10],  This  small  admixture  strongly  affects  the  magnetic 
moment  and  causes  the  ^-factor  of  the  excited  state  to  be  significantly  different  from  that  of  the  ground  state  (~  60%  for  Sr). 
The  different  ^-factor  allows  for  the  addressability  of  the  various  Zeeman  levels  in  the  presence  of  a  bias  magnetic  field  as 
demonstrated  in  Ref.  [25]. 

Optical  spectroscopy  in  atomic  clocks  is  sensitive  to  atomic  motion  due  to  the  Doppler  effect.  To  overcome  this  limitation, 
atoms  are  confined  in  a  tight  optical  lattice  formed  by  a  standing  wave  light  pattern  to  eliminate  broadening  and  frequency  shifts 
due  to  atomic  motion.  Within  a  lattice  site  the  atoms  are  tightly  trapped  in  the  so-called  Lamb-Dicke  regime  where  the  length 
scale  associated  with  their  motion  is  much  smaller  than  the  wavelength  of  the  laser  probing  the  atoms.  Moreover  the  optical 
lattice  can  be  carefully  designed  to  operate  at  the  so-called  magic  wave  length  at  which  the  light  shifts  on  the  clock  states  are 
equal  and  the  clock  frequency  is  not  perturbed  [41]. 

In  optical  lattice  clocks,  one  can  simultaneously  probe  large  samples  ~  103~5  of  laser-cooled  atoms  which  can  potentially 
lead  to  high  frequency  stability.  Nevertheless,  this  precision  may  come  at  the  cost  of  systematic  inaccuracy  due  to  atomic 
interactions.  As  mentioned  before,  the  use  of  ultracold  spin-polarized  fermions  was  thought  to  be  the  key  to  avoid  interaction 
effects.  However,  it  has  been  recently  shown  that  this  is  not  the  case  [11,  13,  19-22,  42], 


B.  Many-body  Hamiltonian  for  spin  polarized  fermionic  atoms 

In  order  to  model  the  dynamics  during  clock  interrogation,  we  will  consider  a  nuclear  spin-polarized  ensemble  of  fermionic 
atoms  with  two  accessible  electronic  degrees  of  freedom  associated  with  the  1Sq-3Po  electronic  levels  (see  Fig.  2)  which  we 
denote  as  g  and  e  respectively,  g  stands  for  the  ground  state  and  e  for  a  excited  state.  We  focus  on  the  case  where  the  atoms 
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FIG.  2.  (Color  online)  Schematics:  a.  Energy  levels  of  alkaline  earth  atoms.  The  1S'o-'iPo  electronic  levels,  which  we  denote  as  g  and  e 
respectively,  define  the  spin-1/2  system.  Here  the  red-blue  colors  indicate  electronic  levels.  The  arrows  show  the  nuclear  spins  which  are  fully 
polarized.  The  indicated  energy  level  spacings  are  specific  of  1 1 1  Yb  atoms.  The  atoms  are  confined  in  a  deep  one-dimensional  optical  lattice 
(generated  by  a  laser  with  wavevector  lattice  along  the  axial  dimension),  which  creates  an  array  of  disk-like  potentials  (pancakes).  In  each 
pancake,  the  atoms  occupy  many  transverse  vibrational  levels  at  current  ~  g K  temperatures,  as  shown  in  c.  The  probing  laser  with  wavevector 
k,  can  be  slightly  misaligned  from  the  lattice  beam,  and  this  gives  rise  to  excitation  inhomogeneities.  The  relevant  interactions  are  p-wave  and 
s-wave.  The  corresponding  interaction  parameters  are  indicated  in  panel  d.  To  derive  the  spin  model,  one  assumes  that  the  atoms  are  frozen  in 
the  vibrational  modes  and  that  the  only  relevant  dynamics,  therefore,  happen  in  the  electronic  degree  of  freedom. 


are  trapped  in  an  external  potential  Vext( R)  which  is  the  same  for  g  and  e  (i.e.  at  the  “magic  wavelength”  [41]).  If  the  atoms 
are  illuminated  by  a  linearly  polarized  laser  beam  with  bare  Rabi  frequency  Q0,  they  are  governed  by  the  following  many-body 
Hamiltonian  [16,  18,  22,  34,  43,  44] 


H  =  H0  +  H1, 

Ho  =  >  Id3 Ri'L (R)  (-£v2  +  ^*(R))  'MR)  +  47r^aeg 
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Here  4fa(R)  is  a  fermionic  field  operator  at  position  R  for  atoms  with  mass  m  in  electronic  state  a  =  g  (1S'o)  or  e  (3Po)  while 
M  R)  =  'H(R)’MR)  is  the  corresponding  density  operator.  We  consider  two  possible  interaction  channels:  s-wave  and 
p- wave  (see  Fig.  2).  Since  polarized  fermions  are  in  a  symmetric  nuclear  state,  their  s-wave  interactions  are  characterized  by 
only  one  scattering  length  a~g,  describing  collisions  between  two  atoms  in  the  antisymmetric  electronic  state,  ^(|<?e)  —  | eg)). 
The  p- wave  interactions  can  have  three  different  scattering  volumes  bgg,  b3e,  and  b3g  associated  to  the  three  possible  electronic 
symmetric  states  (| gg),  |ee),  and  -^=(| ge)  +  \ eg))  respectively.  Hi  takes  into  account  the  interrogation  of  the  atoms  by  a  laser 
that  has  frequency  u>l  and  wavevector  k  and  is  detuned  from  the  atomic  transition  frequency  ojq  by  <5  =  ojl  Co’q. 
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C.  Effective  spin  model 


We  consider  the  situation  of  a  deep  ID  lattice,  Vext  (R),  along  Z,  which  creates  an  array  of  two-dimensional  discs  or  “pan¬ 
cakes”  and  induces  a  weak  harmonic  radial  (transverse)  confinement  with  an  angular  frequency  uir  =  2ttvr.  The  lattice  confines 
the  atoms  to  the  lowest  axial  vibrational  mode.  This  analysis  can  be  straightforwardly  generalized  to  the  case  of  a  2D  lattice  in 
which  two  directions  are  frozen  and  only  one  is  thermally  populated. 

We  expand  the  field  operator  in  a  non-interacting  atom  basis,  \l/a(R)  =  4>q(Z)  can<t>nx  ( X)(j>nY  (F),  where  (f>g  is  the 
ground  longitudinal  mode  in  a  lattice  site  and  <j>n  are  transverse  harmonic  oscillator  eigenmodes.  cj.xn  creates  a  fermion  in  mode 
n  =  ( nx,ny )  and  electronic  state  a.  In  this  basis  and  in  the  rotating  frame  of  the  laser,  H  can  be  rewritten  as  [13,  19,  22]: 

Ho  —  Tie  n  T  ^  *  Enflarn  T  ^  f  (1  ^a,/3  )r£‘S<nn/n/,  n'n  T  V  Hnn'n''nf"  j  danC^n/  C/3n"  Can/,/  5 
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V  Z7T  a 


,3 

a, 
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h ho 
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Here  ?ian  =  4ncan  is  the  atom  number  operator  in  mode  n  and  state  a,  Sa^  is  a  Kronecker  delta  function,  a^0  =  ^Jh/{rnu)} j) 
is  the  radial  harmonic  oscillator  length,  and  En  are  single-particle  energies  in  the  trap.  We  have  used  a  Gaussian  approximation 
for  ()>q  ,  which  is  excellent  for  the  deep  lattice  used  in  experiments,  fnoz  =  2\/ ErVz  with  Er  the  recoil  energy,  H2kfattice/ (2 to), 
^lattice  is  the  lattice  beams’  wave-number  and  Vz  is  the  lattice  depth.  The  coefficients  ,S'rm'n"n"/  and  P, 
and  p-wave  matrix  elements,  respectively,  which  depend  on  the  harmonic  oscillator  modes  and  satisfy  S', 
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Here  Hn(x)  are  Hermite  polynomials. 

In  Fig.  3  we  show  the  mode  dependence  of  the  functions  p(n ,  n',  n,  n')  and  s(n,  n',  n,  n').  Since  those  are  computed  in  the 
harmonic  oscillator  mode  basis,  they  are  long-range  in  mode-space.  While  the  p(n,  n! ,  n,  n')  function  scales  (for  |n  —  n'|  1) 

as  y/n  +  n',  and  grows  with  increasing  energy,  as  expected  from  p- wave  interactions,  the  s(n,n' ,n,n')  function  scales  (for 
| n  —  n'\  1)  as  1/ sj\ n  —  n'\,  and  thus  decreases  with  increasing  energy.  In  Fig.  4,  we  also  show  the  dependence  of  the  mean 

and  standard  deviation  of  the  p-wave  interaction  parameters  Pn n'nn'  as  a  function  of  temperature  (T).  There  one  can  see  that 
Pnn'nn'  is  almost  T  independent  in  this  quasi-2D  geometry.  This  is  expected  because  while  the  actual  p- wave  interactions  for 
fixed  density  should  increase  linearly  with  T  [45],  the  latter  growth  is  compensated  by  the  linear  decrease  with  T  of  the  density 
in  a  2D  harmonic  trap.  Fig.  4  (bottom)  shows  a  histogram  of  Pnn'nn'  which  is  peaked  about  its  average  value.  The  histogram 
was  computed  at  T  =  5/iK  but  it  is  almost  T  independent. 

For  the  laser-atom  interaction  Hamiltonian  we  follow  Refs.  [11,  12,  16-18]  and  assume  that  the  probe  is  slightly  misaligned 
with  a  small  component  along  the  X-direction: 


k  =  k(smipX  +  cos  tpZ), 


(6) 


with  ip  <C  1  the  misalignment  angle  (see  Fig  2).  We  also  assume  we  are  in  a  regime  in  which  laser  induced  sideband  transitions 
can  be  neglected,  and  define  a  mode-dependent  effective  Rabi  frequency  given  by 


fln  —  £loLr 


(v2x)Lo(v2z)e-^+ri-)/2, 


(7) 


jD  ^ 

where  rjx,z  —  kx,zaho  are  Lamb-Dicke  parameters,  a%0  = 


h 

Gjam 


and  Ln  are  Laguerre  polynomials  [46].  In  Fig.  4  we 


show  the  ratio  between  the  mean  Rabi  frequency,  O,  and  the  standard  deviation,  AO,  as  a  function  of  temperature.  AQ  increases 
with  T  since  the  atomic  cloud  spreads  as  it  heats  up.  In  the  Lamb-Dicke  regime,  rjx,z  1,  Hi  becomes: 


H , 


=  51^  4nSen-h.c), 


(8) 


Under  typical  operating  conditions,  vr  ~  450  Hz,  and  as  will  be  discussed  below,  the  system  is  in  the  regime  where  its 
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FIG.  3.  (Color  online)  Mode  dependence  of  the  functions  s(n,  n' ,  n,  n')  and  p(n,  n' ,  n,  n'). 
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FIG.  4.  (Color  online)  The  top  panel  shows  the  mean  (red-dashed  line)  and  standard  deviation  (green-dot-dashed  line)  of  the  p-wave  interactions 
vs  temperature.  It  also  shows  Afl/O  (blue-solid  line)  vs  temperature.  We  used  the  Yb  clock  parameters  with  a  misalignment  of  =  5  mrad. 
The  bottom  panel  is  a  histogram  of  the  Pnn/nn'  at  5  pK.  See  Eq.  (5). 


mean  interaction  energy  per  particle  is  about  two  orders  of  magnitude  weaker  than  the  energy  splitting  between  neighboring 
single -particle  transverse  vibrational  modes.  We  refer  to  this  regime  as  the  vibrational-weakly-interacting  regime.  Thus,  to 
leading  order,  only  collision  events  that  conserve  the  total  single -particle  energy  need  to  be  considered.  Under  these  conditions 
the  atom  population  is  frozen  in  the  initially  populated  modes  and  only  the  electronic  degree  of  freedom  g.  e  vary  during  the 
clock  interrogation.  For  an  initial  state  with  at  most  one  atom  per  mode  (\g) -polarized  state),  it  is  possible  under  these  conditions 
to  reduce  H  to  a  spin-1/2  model  with  the  spin  encoded  in  the  g,  e  states.  For  N  atoms  and  labeling  the  thermally  populated 


(9) 


harmonic  oscillator  modes  as  nj  =  (nxj^Yj)  with  j  G  {1,2,...  iV},  the  spin  model  can  be  written  as: 
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Here  Snj  =  \  Y.a,p  cln. aa/3Cpnr  with  a^,z  Pauli  matrices  and  Inj  =  J2a,p  &ail.cp nj  the  identity  matrix. 
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encapsulate  the  temperature  dependence  of  the  interactions. 

Note,  however,  that  in  the  case  of  a  pure  harmonic  spectmm,  mode  changing  collisions  are  energetically  allowed  even  under 
weak  interactions  due  to  (i)  the  linearity  of  the  harmonic  oscillator  spectmm  and  (ii)  the  separability  of  the  harmonic  oscillator 
potential  along  the  X  and  Y  directions.  Condition  (i)  allows  two  particles  in  modes  (nx  ■  ny)  and  ( mx,nriy )  to  collide  and 
scatter  into  modes  (nx  +  k,  ny  +  k')  and  (rrix  —  k ,  my  —  k')  without  violating  the  energy  conservation  constraint.  Condition  (ii) 
allows  the  same  two  particles  to  scatter  into  modes  (nx,  my)  and  (mx,  fly).  Those  issues,  in  principle,  can  impose  important 
limitations  on  the  validity  of  the  spin  model  in  a  harmonic  trap.  In  practice,  however,  the  trapping  potential  is  not  fully  harmonic. 

—  2  2 

It  comes  from  the  Gaussian  beam  profile  of  the  lasers  and  is  given  by  Vr  ~  —Ae  with  vjq  the  beam  waist.  To  leading  order, 

2 

the  trapping  potential  is  harmonic  Vr  ~  R2  but  for  an  atom  in  a  mode  {nx,  ny}  there  are  higher  order  corrections  of  the 
energy  beyond  leading  order:  En  =  fkOR(nx  +  ny  + 1)  +  A En,  with  A En  ~  Hlor,  ^3( n2x  +»y)+  4 nxny  +5 (nx  + 

ny  +  1)^ .  At  typical  operating  conditions  of  the  Yb  and  Sr  experiments:  uir  ~  27r  x  450Hz,  Wq  ~  30  —  100/im,  T  >  1/iK, 
and  a  mean  occupation  mode  numbers  nx.y  >  50,  the  difference  of  A En  for  nearby  modes  is  larger  than  27 r  x  10  Hz  which  is 
not  negligible  compared  to  typical  interaction  energy  scales  -'-Hz.  The  first  term  in  A En  thus  prevents  processes  (i),  while  the 
second  term,  which  breaks  the  separability  of  the  potential,  prevents  processes  (ii).  Based  on  this  argument  we  first  restrict  our 
analysis  to  only  processes  that  conserve  the  number  of  particles  per  mode. 

A  further  simplification  of  Eq.  (9)  can  be  made  when  atoms  are  initially  prepared  in  the  totally  symmetric  Dicke  manifold 
with  S  =  N/2  [47]  at  time  t  =  0.  Here,  S(S  +  1)  is  the  eigenvalue  of  the  collective  operator  S  ■  S  and  5,T=a 
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In  this  case  there  are  two  important  physical  mechanisms  that  prevent  leakage  of  the  population  outside  the  symmetric 
Dicke  manifold,  (i)  The  weak  dependence  of  the  interaction  matrix  elements  on  the  thermally  populated  modes  so  that  the 

Am  "  ■ 


mode-dependent  coupling  constants  n,  ,Stlj ,  \rl  ,  and  Cn  are  peaked  at  their  averages  JT  = 


N(N-l) 


Xn 


5 C'n 


— —  and  Ca  =  (See  Figs.  3-4).  (ii)  The  fact  that  jA  »  Axh,  ACg.  Here  Axh,  ACs  are  the  corre¬ 

sponding  standard  deviations.  The  latter  is  satisfied  in  part  because  JA  is  the  only  interaction  term  that  has  a  contribution  arising 
from  s-wave  interactions.  In  general,  those  are  expected  to  dominate  over  p- wave  interactions  since  p- wave  collisions  are  sup¬ 
pressed  by  the  centrifugal  barrier  which  is  estimated  to  be  greater  than  ~  25/iK  [11,  19].  It  must  be  said,  nevertheless,  that  the 
actual  values  of  the  s-wave  and  p-wave  scattering  parameters  between  two  e  atoms  and  one  e  and  one  g  are  not  known,  (i)  and 
(ii)  impose  a  large  energy  gap  in  the  Hamiltonian  which  suppresses  transitions  between  manifolds  with  different  total  collective 
spin  S,  caused  by  the  inhomogeneities  A\'  and  AC  [48].  Consequently,  to  a  very  good  approximation  it  is  expected  that  the 


FIG.  5.  (Color  online)  In  Ramsey  spectroscopy  two  short  pulses  of  area  6\  and  do  are  applied,  separated  by  a  dark  time  r. 


dynamics  can  be  projected  into  the  S  =  N/2  manifold  with  an  effective  Hamiltonian  given  in  terms  of  collective  operators: 

H$  =  HS  +  Hg  (16) 

Hs/h  =  - 6SZ  +  J^S-S  +  xa  (SZY  +  Ca(N-  1)  Sz, 

H$  =  - nnsy ,  (17) 

The  term  S  ■  S  is  a  constant  of  motion  and  does  not  play  any  role  in  the  collective  dynamics. 

In  addition  to  the  above  conditions  on  the  interactions,  staying  in  the  S  =  N/2  manifold  requires  that  the  probing  laser 
generates  negligible  excitation  inhomogeneity,  i.e.  A 0/0  <c  1.  In  this  case  we  can  write  the  atom-light  Hamiltonian  in  terms  of 
collective  operators.  The  validity  of  this  condition  depends  on  the  misalignment  angle,  tp  [Eq.  (6)],  and  the  average  vibrational 
mode  population  [12],  At  the  current  operating  conditions,  the  JILA  Sr  clock  can  achieve  AO/O  <0.1  [22],  and  the  collective 
mode  approximation  for  the  atom-light  Hamiltonian  is  well  satisfied.  The  latter  is  not  necessarily  the  case  for  the  Yb  optical 
clock  experiment  at  NIST  due  to  the  higher  temperatures,  as  we  will  elaborate  later.  The  restriction  of  the  dynamics  to  the  Dicke 
manifold  is  relaxed  in  Sec.  VI  where  we  also  investigate  the  parameter  regime  in  which  it  is  valid. 


III.  RAMSEY  INTERROGATION:  COLLECTIVE  CASE 


In  Ramsey  spectroscopy  (see  Fig.  5),  a  well  established  tool  in  atomic  physics,  atoms  are  typically  prepared  in  the  same 
internal  state,  say  g  (e.g.  via  optical  pumping).  Next,  one  applies  a  strong  resonant  linearly  polarized  light  pulse  for  time  t\. 
By  strong,  we  mean  that  the  Rabi  frequency  must  be  much  larger  than  the  atomic  interaction  energy  scales  but  weaker  than  the 
harmonic  oscillator  frequency,  vr,  to  avoid  laser  induced  mode  changing  processes. 

This  first  pulse  rotates  the  spin  state  of  the  atom  at  mode  n  to  cos(0")|g)  +  sin(6l“) | e),  here  is  the  pulse  area. 

Subsequently,  the  atoms  are  allowed  to  freely  evolve  for  a  dark  time  r.  Finally  a  second  pulse  is  applied  for  a  time  £2  and  the 
population  of  the  e,  g  states  measured. 

In  the  following,  we  will  discuss  the  case  of  a  fixed  number  of  atoms,  N,  prepared  in  the  modes  n  =  {ni , . . . ,  iin},  and  only 
later  we  will  worry  about  spatial  and  thermal  averages.  Note  that  under  the  collective-mode  approximation  9 1,2  =  fVi  .2- 


A.  Analytic  solution 


Let  us  first  start  dealing  with  the  situation  in  which  there  is  no  excitation  inhomogeneity  and  the  initial  (second)  pulses  are 
just  rotations  of  the  collective  Bloch  vector  by  an  angle  9\  (9 2),  respectively.  In  this  case,  the  state  of  the  system  before  the  free 
dynamics  is 


wn>  =  E 


o\ 


(cos—) 


c(sin^)fc| N/2, k- N/2), 


with  | N/2,  M )  being  collective  Dicke  states.  During  the  free  evolution,  the  Hamiltonian  reduces  to  Xa(Sz)2  +  Ca(N  —  1  )SZ 
which  introduces  just  a  phase,  e~lT(XR(M'>2+Cs('N^r>M^  to  each  of  the  | N/2,  M)  states.  Here  M  is  an  eigenvalue  of  Sz  and 
takes  integer  or  half-integer  values  (depending  on  whether  N  is  even  or  odd)  satisfying  — N  <  2 M  <  N.  Expressions  for  the 
evolution  of  the  spin  operators  after  a  dark-time  evolution  r  can  be  exactly  computed: 
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with  Z„  and  £a  and  A given  by 
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{Sy) 
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N  cos  0i 
2  ’ 


(18) 

(19) 

(20) 


Z?  =  1  -  sin2  (0i )  sin2  (xa t)  , 

(21) 

tan(Cn-r)  =  tan  (xht)  cos(0i), 

(22) 

27rAzzfi  =  (N  -  l)(C7a  -  Ca). 

(23) 

The  normalized  contrast,  which  is  a  measure  of  the  amplitude  of  the  Ramsey  fringes,  is  defined  as  the  magnitude  of  the 
projection  of  the  collective  Bloch  vector  on  the  x-y  plane  normalized  by  half  of  the  total  number  of  atoms: 


from  the  above  expression  one  obtains  that 


C  = 


2\J  (Sx)2  +  (Sy)2 
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C  = 


sin  (0JZ*-1 


(24) 


(25) 


The  contrast  is  extracted  from  measurements  of  the  fraction  of  excited  atoms,  iVe,n(ti,  (2) /N  by  varying  the  laser  detunning, 
or  by  varying  the  phase  of  the  second  pulse  along  the  x-y  plane.  For  the  former  case: 


C  = 
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Ne. 


+  Ne. 


-  2 (Ve,fi 

t5=7V  tS  =  7T /2/ 
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(27) 


Let’s  now  discuss  the  physics  encapsulated  in  Eqs.  (18-20).  The  term  Av^  is  the  so  called  density-dependent  frequency  shift,  which  gives 
rise  to  a  density-dependent  measurement  of  the  atomic  transition  frequency.  The  quantity  Z^  ~ 1  determines  the  contrast  of  the  Ramsey  fringes. 
In  the  weakly  interacting  regime,  xrt  "C  1,  and  for  N  1,  27rAzzg  «  ( NCa  +  2xn(>S2}).  This  means  that  interactions  act  as  an  effective 
magnetic  field  along  the  quantization  axis  with  magnitude  £>eff  =  (NCa  +  2xa(Sz)),  which  depends  both  on  the  total  atom  number  and 
the  population  difference  between  excited  and  ground  state.  This  is  consistent  with  just  a  simple  interpretation  of  the  frequency  shift  as  being 
the  average  energy  difference  experienced  by  an  atom  in  state  e  with  respect  to  an  atom  in  state  g  due  to  the  presence  of  other  atoms.  As  we 
explain  in  Sec.  Ill  B,  this  can  be  derived  from  a  mean  field  analysis  which  factorizes  the  interaction  term  as  {(Sz)2)  oc  SZ(SZ).  In  this  regime, 
the  condition  cos(0°pt)  =  Ca/x fl  determines  the  pulse  area  0°pt  at  which  the  shift  is  canceled.  This  is  the  ideal  operating  pulse  area  for  a 
clock  [20],  Note  that  if  Ca  is  equal  to  zero,  no  density  shift  is  expected  at  ( Sz )  =  0.  This  is  consistent  with  the  intuition  that  due  to  the  equal 
population  of  both  e  and  g  states,  the  mean  energy  shift  experienced  by  an  atom  in  state  e  due  to  others  is  exactly  canceled  by  the  opposite 
energy  shift  experienced  by  an  atom  in  g.  In  this  weakly  interacting  regime,  Za  ~  1  as  expected  in  the  case  that  interactions  act  just  as  a  mere 
effective  magnetic  field,  which  causes  the  Bloch  vector  just  to  precess  with  no  Ramsey  fringe-contrast  decay. 

Outside  the  weakly  interacting  regime,  two  important  corrections  to  this  picture  arise.  One  is  the  fact  that  the  shift  is  no  longer  linear  in  ( Sz ) 
[from  Eq.  (22)].  The  second  one  is  that  the  Ramsey  fringe-contrast  collapses  and  revives.  The  collapse  is  well  approximated  by  a  Gaussian 
decay,  Z A  ~  e~N/ 2sm  (9l)*aT  .  The  revivals  take  place  at  times  XnT„  =  nn  with  n  an  integer.  This  behavior  of  the  contrast  is  closely 
linked  to  the  decay  of  coherence  in  a  matter-wave  due  to  the  nonlinearities  arising  from  the  atom-atom  interactions  and  subsequent  revival  due 
to  the  discreteness  of  the  spectrum  of  the  many-body  system.  The  observation  of  collapses  and  revivals  of  a  Bose-Einstein  condensate  (BEC) 
loaded  in  an  optical  lattice  was  first  reported  in  Ref.  [49]. 


B.  Mean-field  solution 

Even  though  the  all-to-all  interactions  allows  for  an  exact  solution  of  the  many-body  dynamics  in  Ramsey  spectroscopy,  it  is  convenient 
to  introduce  an  approximate  mean-field  treatment.  The  mean-field  treatment  will  be  very  helpful  for  dealing  with  inelastic  collisions,  which 
are  experimentally  relevant.  A  simple  and  enlightening  way  to  carry  out  a  mean-field  treatment  is  to  use  the  Schwinger-boson  representation, 
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FIG.  6.  (Color  online)  Comparisons  between  the  GPE  (mean-field,  red-dashed  line),  the  TWA  (  green-dot-dashed  line)  and  the  exact  solution 
(black-solid  line).  Here  we  assumed  N  =  30  and  Cr  =  5  =  0.  The  top  plot  is  for  6 1  =  7t/2,  the  lower  plot  for  0\  =  7t/4. 


which  maps  spin  operators  to  two-mode  bosons  subject  to  a  constraint  [50],  It  represents  the  spin  operators  as 

2  sz  =  4'1'I'e  - 

S+mVt&g, 

§-  = 


(28) 

(29) 

(30) 

(31) 


with  the  constraint 


N  = 


(32) 


with  (E'a  a  bosonic  annihilation  operator  of  mode  a  =  e,g.  The  mean-field  treatment,  which  gives  rise  to  the  so-called  Gross-Pitaevskii 
Equation  (GPE)  [51],  replaces  the  field  operators  Tc.g  by  c-numbers,  'Pe,9  —¥  The  latter  approximation  is  justified  in  the  weakly 

interacting  regime  when  there  is  a  macroscopic  population  of  those  modes  and  the  state  of  the  system  can  be  regarded  as  a  simple  product 
state  with  no  significant  entanglement.  Translated  to  the  spin  language,  those  conditions  imply  that  the  system  can  be  well  described  as  a  spin 
coherent  state.  For  purposes  that  will  become  clearer  later,  it  is  convenient  to  introduce  the  density  matrix  paa  =  Na  and  pap  =  The 

latter  satisfies  the  following  equations  of  motion: 


^(5z)gpe  =  ti(Sx)GPE, 


d  d 

Pgt 

Meg  =  -ti{Sz)GFE  -  i  1(5  -  NCr  -  2Xn(<§z)GPE| 


dtPe9  ~  d/3e  ~  Me9’ 


peg- 


(33) 

(34) 


(35) 
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Note  pge  =  plg  and 


{SZ)GPE 


(SX)GPE 


(SV)GPE 


pee  pgg 
2 

peg  +  Pge 
2 

Peg  —  Pge 

2  i 


(36) 

(37) 

(38) 


are  the  components  of  the  Bloch  vector. 

For  the  Ramsey  dynamics,  the  mean-field  treatment  is  almost  trivial,  and  time  evolution  corresponds  to  a  precession  of  the  Bloch  vector 
induced  by  an  effective  magnetic  field 

Bf  =  2nAu§PE  =  V[Cfl  -  Xfl  caKflr)].  (39) 


This  behavior  can  be  clearly  seen  in  Eq.  (34)  by  noticing  that  during  the  dark  time  fi  =  0,  thus  (S'Z)GPE  is  a  constant  of  motion,  and 
peg  precesses  at  a  rate  8  —  2-kAv<~pe‘.  Note  that  in  striking  disagreement  to  the  exact  many-body  solution,  at  mean-field  level,  the  Ramsey 
fringe  contrast  never  decays  [52].  This  implies  that  it  agrees  with  the  exact  solution  only  when  Z ~  -  1  <  1,  which  can  be  satisfied  when 
y'N  sin(#i)xni"  1  ,  i.e.  short  times,  weak  interactions,  or  small  pulse  areas. 


C.  Truncated  Wigner  Approximation 


The  truncated  Wigner  Approximation  (TWA)  has  proven  to  be  a  successful  approach  to  incorporate  the  leading  quantum  corrections  to  the 
mean-field  dynamics.  See,  for  example,  Refs.  [26,  27]  and  references  therein.  To  implement  the  TWA,  one  needs  to  solve  the  mean-field 
equations  of  motion  supplemented  by  random  initial  conditions  distributed  according  to  the  Wigner  function.  For  a  spin  coherent  state  with 
S  =  N/ 2,  pointing  initially  in  the  x-z  plane,  hg1  =  (sin  9i ,  0,  —  cos  9 1 ) ,  the  Wigner  function  is  given  by: 


p(So,Sq,Sq)  =  (  —  )  <5(So  cos 9\  -  Sg  sin(9i  +  S)e 


(Sq  sin  Ai+Sq  ■ 


3*ir 


(40) 


This  Wigner  function  has  a  transparent  interpretation.  For  example,  if  the  Bloch  vector  points  along  the  direction  —z,  i.e.  9 1  =  0,  because  of 
the  uncertainty  principle,  the  transverse  components  still  fluctuate  so  that  (Sv Sv)  =  { SXSX )  =  S/2.  Quantum  mechanical  expectation  values 
of  the  spin  operators  relevant  for  this  work  can  be  computed  as  ( 0(t ))  =  f  dSgdSgdSg  ( 0(r))p(Sg ,  Sg,  Sg)  with  (0(t))  the  classical 
evolution  of  the  observable  calculated  using  the  mean-field  equations. 

Using  this  prescription  we  obtain: 


(Sx(t))  =  Re[(S+(r))], 

(41) 

{Sy  (t)>  =  Im[(S+(r))], 

(42) 

(43) 

<S+(t))  =  Arsin(^1)(i  +  iXf.7.cos[01])e-§x«T2sin2(0l)eiT(5-27rA‘/*PE). 

(44) 

Flere  Re  and  Im  correspond  to  the  real  and  imaginary  parts,  respectivelly.  The  essential  physics  is  encapsulated  in  the  term  that  exhibits  a 
Gaussian  decay,  proportional  to  e-TXfiT  sln  and  in  the  phase  shift,  proportional  to  2nAi/?PE .  The  exponential  term  leads  to  a  decay  of 
the  Ramsey  fringe  contrast  missing  at  the  mean-field  level.  Note  that  in  addition  to  the  exponential  decay  of  the  contrast  and  the  mean-field 
phase  shift,  there  is  an  extra  term,  ixfii" cos[0i]  in  Eq.(44).  For  |XnTl  <  1,  this  term  just  introduces  a  finite  N  correction  (i.e.  replaces 
N  — >  N  —  1)  in  the  phase  shift,  irrelevant  in  the  large- At  limit.  The  TWA  approximation,  nevertheless,  fails  to  reproduce  the  periodic 
revivals  of  the  coherence,  which  take  place  in  the  exact  solution.  Ftowever,  for  current  experimental  probing  times,  which  are  much  shorter 
than  the  revival  time,  the  TWA  is  an  excellent  alternative  for  capturing  the  quantum  dynamics.  Fig.  6  shows  comparisons  between  the  GPE, 
the  TWA  and  the  exact  solution  for  different  pulse  areas  and  for  N  =  30.  The  power  of  the  TWA  is  demonstrated  in  its  capability  to  fully 
capture  the  decay  of  the  Ramsey  fringe  contrast.  We  note  that  other  approximation  methods,  typically  used  to  account  for  quantum  fluctuations 
beyond  mean-field,  such  as  the  time-dependent  Bogoliubov  approximation,  generally  fail  to  reproduce  the  decay  of  quantum  coherences  [53]. 
Moreover,  as  we  will  show  below  in  Sec.  IV,  the  TWA  can  be  straightforwardly  generalized  to  deal  with  the  many-body  dynamics  of  an  open 
quantum  system. 
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FIG.  7.  (Color  online)  Comparisons  between  the  GPE  (red-dashed  line),  the  TWA  (green-dotted  line)  and  exact  solution  (black-solid  line). 
Here  we  assumed  N  =  30,  C„  =  5  =  0,  Fg’e  =  Xfl/3,  and  Ft’9  =  0.  The  plots  are  for  9i  =  ir /2  (left),  for  9\  =  7t/4  (middle)  and 
9i  =  37r/4  (right). 


IV.  RAMSEY  INTERROGATION:  INCLUDING  LOSSES 


A.  Master  equation 


The  Hamiltonian  formulation  described  above  is  valid  only  for  a  closed  system.  To  account  for  losses  due  to  inelastic  e-e  or  e-g  collisions, 
recently  observed  in  experiments  [20,  21],  one  needs  to  use  instead  a  master  equation: 


h^p  =  -i  [Hi,p\  +  Cp. 


(45) 


Here  p  is  the  reduced  density  matrix  operator  of  the  many-body  system.  is  the  Hamiltonian  given  by  Eq.  (9),  and  £  is  a  is  a  Lindbladian 
superoperator  that  accounts  for  inelastic  processes.  Considering  p-wave  e-e  and  e-g  losses  (s-wave  losses  do  not  occur  if  one  restricts  the 
dynamics  to  the  fully  symmetric  manifold  [54]  as  we  will  assume  below),  C  is  given  by 


C 

h 


5Z  rib’IV  )+] 

j^j'  ,a.=ee,eg 


j^j'  ,cx=ee,eg 


)tAnj,n,,P  +  p(A n  , 


(46) 


Here  the  jump  operators  are  y4„e  n  t  = 
expression  for  7ee’eff  is  identical  to  ve,e 
the  inelastic  one. 


-  Cen.Cen.,  and  A^3  =  (cen,cgn.,  +  cgn . cen  , )/V2,  while  =  -yee.e9pn ;n ., ■  The 

•J  J  3  '  J  4  J  J  J  J  5  J  J  J 

and  ve,a  respectively  (see  Eq.  (3))  up  to  the  replacement  of  the  p-wave  elastic  scattering  volume  by 


The  last  two  terms  in  Eq.  (46)  can  be  thought  of  as  generating  an  effective  Hamiltonian,  in  which  the  terms  {V^f  n ,  V"n9.n  , }  are  replaced 
by  {(I'ri'jn  /  —  frne,n.,)5  {Vn3,n  ,  ~  ,)}  respectively.  On  the  other  hand,  the  first  term  in  Eq.  (46),  the  so-called  recycling  term, 

requires  a  full  density  matrix  formulation  and  thus  the  dimension  of  the  Hilbert  space  grows  much  faster  with  N  than  the  corresponding 
dimension  of  the  pure  Hamiltonian  case.  For  example  just  for  4  particles  the  required  Hilbert  space  is  ~  103  states.  Nevertheless,  keeping  the 
first  term  is  crucial  to  conserve  the  trace  of  the  density  matrix. 

An  important  point  to  keep  in  mind  is  that  the  recycling  term  connects  sectors  of  the  Hilbert  space  with  different  atom  number.  However, 
since  the  environment  knows  how  many  particles  left  the  system,  there  are  no  coherences  between  sectors  of  different  atom  numbers  and  the 
master  equation  can  be  solved  in  a  “block-diagonal  way”  [55], 

The  master  equation  can  be  further  simplified  under  the  collective  mode  approximation.  Let  pjj  be  the  density  matrix  for  a  single  sector  of 
N  particles,  and  let’s  assume  that,  at  time  t  =  0,  we  start  in  the  sector  with  N  =  N  particles.  As  explicitly  shown  in  Appendix  1,  to  solve 
for  the  master-equation  dynamics,  we  need  to  solve  a  series  of  differential  equations  for  each  of  the  subspaces  with  cascading  atom  numbers. 
The  sector  Af  =  N  does  not  have  driving  terms  and  can  be  solved  by  merely  solving  the  effective  Hamiltonian  dynamics.  We  then  use  the 
integrated  solutions  in  the  A f  =  N  sector  as  a  driving  terms  for  Af  =  N  —  2,  the  latter  for  Af  =  N  —  4,  etc.  After  the  dynamics  of  all  the 
density  matrix  sectors  are  known,  one  can  compute  any  observable.  For  N  <  50,  the  above  procedure  can  be  efficiently  performed  numerically 
|56|. 
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B.  Mean-field  Treatment 

At  the  mean-field  level,  accounting  for  losses  is  significantly  simpler.  We  assume  that  the  reduced  density  matrix  of  the  many-body  system 
can  be  factorized  as: 

/5  =  0pflj  pa3 ■=  Pa,p\a)W  (47) 

j  cx,f3=e,g  ,0 

Here  pa.  is  the  reduced  density  matrix  of  atom  in  mode  n, .  Here  g ,  e  label  the  two  possible  spin  states  of  the  atom  and  0  is  the  vacuum. 

By  using  this  ansatz  in  Eq.  (46),  one  can  write  closed  equations  of  motion  for  (Pa  p.  At  the  mean-field  level,  the  equations  of  motion  for 
Pa  p  f°r  a>  P  £  (e,  <?)  are  fully  decoupled  from  the  ones  for  p3a  0  and  pg  a,  meaning  that  if  there  are  not  particle-hole  coherences  initially  in 
the  system,  they  do  not  develop  in  the  mean-field  dynamics.  Thus  the  only  role  the  recycling  terms  play  in  the  dynamics  is  to  populate  the 
vacuum  state,  which  can  be  accounted  for  by  the  constraint:  Yla=e  g  Pa, <x  =  1  —  Po  o- 

When  the  recycling  terms  do  not  contribute  to  the  equations  of  motion  of  pJa  /3.  the  effect  of  the  losses  can  be  incorporated  by  adding  to 
the  Hamiltonian  a  non-hermitian  effective  part.  For  the  case  of  all-to-all  interactions,  the  total  net  Hamiltonian  that  determines  the  mean- 
field  dynamics,  including  single-particle  terms  and  both  elastic  and  inelastic  collisions,  can  be  written  in  terms  of  the  collective  variables 

EN 

j  =  1  Pa,p  as: 


HaF/H  =  7^  (Pee  Pgg )  -  ^  (peg  ~  Pge)  ~  ^f{Pee  ~  Pggf  +  ^yOee  ~Pgg) 

pe,e  pe,g 

-i-j-Pee  ~  i-^-peePgg, 


with  rj®  = 


and  Ti’e  = 


N(N-l)  r,  ~  N(N-l) 

The  collective  density  matrix  equations  of  motion  obtained  from  that  Hamiltonian,  Eq  (48)  are: 


(48) 


d  _  H  x  pe.g 

Qf  Pg3  —  2  ^e®  +  ^®e  PeeP gg> 

JJJ  P<i‘  =  ~^{peg  +  Pge)  -  pee  —  r  3®PeePgg, 


d 

-Q^  PeS  ~~  Mes  ~ 


rr^e  ,  ^  "a  {Pee  +  Pgg) 


peg- 


(49) 

(50) 

(51) 


Note  the  term  Meg  gives  the  mean-field  dynamics  of  peg  in  the  absence  of  losses  and  was  defined  in  Eq.  (35).  The  terms  proportional  to  Ti’® 
and  rt’ e  give  rise  to  classical  rate  equations  for  the  population  in  the  presence  of  two-body  losses. 


C.  Truncated  Wigner  Approximation 

Following  a  similar  procedure  to  the  one  explained  for  the  Hamiltonian  case  without  decay,  one  can  include  quantum  correlations  by  first 
solving  the  mean-field  equations  of  motion  determined  by  Eqs.  (49-51)  supplemented  by  random  initial  conditions  distributed  according  to 
the  Wigner  distribution,  Eq.  (40).  In  Fig.  7,  we  compare  the  solution  of  the  exact  many-body  master  equation,  the  mean-field,  and  the  TWA 
with  losses.  We  can  see  some  difference  between  the  TWA  and  the  master  equation  for  N  =  30,  but  overall  the  TWA  does  an  excellent  job 
reproducing  the  full  dynamics. 


V.  IMPROVED  SPIN  MODEL:  INCLUDING  VIRTUAL  MOTIONAL  EXCITATIONS 

A.  Resonant  terms 


As  discussed  in  Sec.  II C,  in  a  pure  harmonic  spectrum,  mode  changing  collisions  are  energetically  allowed  and  impose  important  limitations 
on  the  validity  of  the  spin  model.  However,  when  one  takes  into  account  the  Gaussian  shape  of  the  actual  potential,  corrections  to  the  harmonic 
spectrum  are  at  the  level  of  the  interaction  energy  at  current  experimental  conditions  and  are  enough  to  prevent  mode-changing  collisions. 
Consequently,  resonant  transitions  arising  from  the  perfect  linearity  of  an  harmonic  oscillator  spectrum  can  be  treated  as  off-resonant.  Off- 
resonant  terms,  nevertheless,  can  provide  corrections  to  the  dynamics  predicted  by  the  spin  model,  and,  in  this  section,  we  explain  a  way  to 
incoiporate  those  in  the  dynamics. 
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B.  Off-Resonant  terms 

The  spin  Hamiltonian  neglects  collision  processes  that  do  not  preserve  the  single-particle  energy.  However,  off-resonant  collisions  can 
still  take  place  virtually  and  will  introduce  corrections  to  the  spin  model.  To  account  for  those,  we  split  the  N-particle  Hilbert  space  into  the 
resonant,  E,  and  off-resonant,  T,  manifolds  respectively,  spanned  by  the  states: 

N 

'nj,  (52) 

3  =  1 
N 

Etot  =  J2E^ +Kot- 

2=1 

These  states  are  written  in  the  occupation  basis  and  a  £  { g ,  e}.  Since  the  same  mode  kj  can  be  occupied  simultaneously  by  a  g  atom  and  an 
e  atom,  gk1  and  ek,  can  occur  simultaneously  in  T  J  ) . 

3  J  k 7 

The  spin  Hamiltonian  was  obtained  by  directly  projecting  the  interaction  part  of  the  many-body  Hamiltonian  on  E.  An  effective  Hamiltonian 
that  accounts  for  the  leading  order  corrections  generated  by  virtual  occupation  of  off-resonant  states  can  be  derived  by  using  the  Schrieffer- 
Wolff  transformation  [57]: 

Hf  =  Hi  +  H ?2.  (53) 

Here,  H§  is  the  spin  model  given  by  Eq.  (9).  FT?2  is  obtained  via  second-order  perturbation  theory  as  follows: 


—  I^ni  5  crn2  5  •  • 

•  ?  ^"n  tv 

vpT  ) 

—  I^ki  5  ^k2  j  •  • 

y,  „  Q  „  (<I>5 jH\*l .)<*£ »  ) 


tot  _ 

k  0 


(54) 


Here  H  is  the  interaction  part  of  Eq.  (3). 

Since  if?2  acts  on  the  E  subspace,  it  is  convenient  to  explicitly  write  it  in  terms  of  spin  operators.  To  accomplish  that,  we  divide  the  states 
|4>J.)  into  6  different  categories,  |\Kl£),  ),  |T'i'[ld),  which  are  shown  in  Fig.  8. 

k  fc  k  k  k  k  k 

The  states  labeled  with  superscripts  1  h  and  2h2m  are  states  which  have  one  and  two  atoms  populating  a  mode  not  belonging  to  the  initially 
populated  n  =  {ni,  112, . . . ,  njv}  manifold  respectively.  The  states  labeled  with  superscripts  2hlm  are  states  which  have  two  atoms  both 
populating  one  mode  not  belonging  to  the  initially  populated  manifold.  The  states  labeled  with  superscripts  Id  and  2d  have  one  and  two  modes 
in  the  initially  populated  manifold  occupied  by  two  atoms  (doublon)  respectively.  Finally  states  labeled  by  superscripts  lhld  have  one  atom 
populating  a  mode  outside  the  initially  occupied  manifold  and  a  doublon.  We  will  consider  the  various  contributions  case  by  case  in  Appendix 
2.  Here  we  only  quote  the  final  result,  assuming  collective  interactions  and  dropping  terms  that  are  constants  of  motion. 


FIG.  8.  (Color  Online)  Schematic  representation  of  the  various  virtual  processes.  Here  the  filled  circles  denote  g  (blue)  and  e  (red)  atoms  in 
the  initially  populated  levels  (also  shown  with  a  solid  line).  The  striped  circles  denote  atoms  virtually  populating  off-resonant  levels,  shown 
with  a  dashed  line.  Note  that  the  classification  is  made  in  terms  of  the  number  of  doublons  or  populated  states  outside  the  initially  populated 
manifold.  Processes  in  which,  in  addition  to  the  formation  of  a  hole  or  a  doublon  within  the  occupied  manifold,  the  spin  is  switched  between 
two  atoms  are  also  included  in  each  classification  but  are  not  necessary  shown  in  the  schematics.  Also,  since  this  is  a  schematic  plot,  for 
simplicity,  we  show  two  e  atoms  and  two  g  atoms  only  and  the  color  was  selected  randomly.  For  example,  for  the  cases  (lh),  (Id),  (2d)  and 
(lhld)  corresponding  processes  exist  where  g  (blue)  atoms  are  the  ones  that  get  excited. 
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M2  =alSz +al{Szf +  al{Szf,  (55) 

The  parameters  a^2,3  are  explicitly  computed  in  Appendix  2. 

Among  those  corrections,  linear  and  quadratic  contributions  can  be  absorbed  into  the  spin  Hamiltonian,  and  only  the  processes  proportional 
to  (52)3  give  rise  to  an  additional  term.  Multi-body  interactions  arising  from  virtual  scattering  processes  have  been  shown  to  introduce 
measurable  corrections  in  quantum  phase  revivals  in  the  Bose-Hubbard  model  [58,  59],  The  cubic  terms  can  give  rise  to  visible  corrections  to 
the  Ramsey-fringe  contrast  as  recently  measured  in  the  JILA  Sr  clock  and  reported  in  Ref.  [22],  Here  we  also  show,  in  Sec. VII B  2,  that  they 
play  an  important  role  in  the  Yb  Ramsey-fringe  contrast  decay. 

In  principle,  terms  proportional  to  (S'2)4  could  appear  in  the  effective  Hamiltonian.  We  find,  however,  that  they  vanish  in  our  system  due 
to  the  fact  that  when  an  atom  in  mode  tij  scatters  to  an  already  occupied  mode,  n , / ,  forming  a  doublon,  the  reverse  process  (an  atom  in  mode 
n j>  scatters  to  mode  n,  forming  a  doublon),  contributes  with  the  opposite  sign  to  the  effective  model,  and  thus,  both  cancel. 


VI.  RAMSEY  INTERROGATION:  BEYOND  THE  COLLECTIVE  REGIME 

Negligible  excitation  inhomogeneity  is  not  necessarily  a  good  approximation  in  all  optical  lattice  clock  experiments.  For  example,  in  the 
Yb  optical  lattice  clock  at  NIST,  the  current  operating  temperature  is  typically  around  10/iA'  and  A Q/Cl  can  be  in  some  cases  as  high  a  0.25 
(see  Fig.  4).  In  the  following  sections  we  describe  various  ways  to  treat  inhomogenous  excitation  and  interaction  effects  that  cause  population 
of  non-collective  states.  We  first  consider  the  case  of  collective  interactions,  where  an  analytic  solution  to  the  full  quantum  problem  in  the 
presence  of  an  excitation  inhomogeneity  can  be  derived  based  on  perturbation  theory  in  Afi/fi.  Next,  we  discuss  a  solution  based  on  the  TWA 
that  incorporates  in  the  model,  not  only  excitation  inhomogeneity  but  also  inhomogeneity  in  the  spin  coupling  constants  and  losses. 


A.  Analytic  solution 

In  the  presence  of  excitation  inhomogeneity,  an  analytic  treatment  based  on  perturbation  theory  can  be  performed  under  the  collective 
Hamiltonian  approximation,  i.e.  neglecting  two  body  losses  and  assuming  collective  two-body  interactions. 

To  accomplish  that  one  writes  finj  =  fig  +  ,  with  Dg  =  JY  flnj  /N  the  mean  Rabi  frequency  and  treats  SQ.nj  /Cl  as  a  perturbation 

parameter.  Keeping  only  the  leading  order  correction,  which  can  be  shown  to  be  quadratic  in  AS2n  =  ^^Y  Gnj/IV  —  i.e.  the  first-order 
corrections  vanish,  and  after  a  lengthly  but  straightforward  calculation  described  in  Appendix  3,  one  can  compute  the  number  of  excited  atoms, 
£2),  in  Ramsey  spectroscopy.  The  expression  is  quite  complicated  and  we  do  not  explicitly  show  it  in  the  main  text.  Instead  here  we 
discuss  its  behavior  in  the  weakly  interacting  regime.  \rt  <C  1  and  1,  where  one  can  show  that  the  density  shift  becomes: 
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with  A 0"2  =  Af)flfi,2  and  0]  .2  = 

The  expression  of  the  shift  in  the  weakly  interacting  limit  is  very  illuminating.  Note  that  while  to  leading  order  (AO  —  0)  the  density  shift  is 
only  determined  by  the  excitation  fraction,  cos(Rl)  =  (Ng  —  Ne)/ 2,  and  is  independent  of  the  second  pulse  area,  8  2.  excitation  inhomogeneity 
does  introduce  a  dependence  of  the  density  shift  on  the  second  pulse  area.  An  intuitive  explanation  for  that  can  be  obtained  by  noticing  that  the 
density  shift  is  extracted  from  a  collective  measurement.  Therefore,  if  there  are  atoms  excited  outside  the  symmetric  Dicke  manifold  due  to 
the  first  pulse,  they  only  contribute  to  the  measured  signal  if  they  are  brought  back  to  the  symmetric  Dicke  manifold  by  the  second  pulse.  The 
mean-field  expression,  Eq.  (56),  also  tells  us  that  terms  proportional  to  (AO" )2  just  slightly  renormalize  the  overall  magnitude  of  the  density 
shift,  Xfl  — >  Xit’  without  affecting  its  general  dependence  on  the  excitation  fraction.  On  the  other  hand,  the  term  proportional  to  AO"  AO" 
does  modify  the  dependence  of  the  density  shift  on  pulse  area  and  vanishes  at  0J  =  7r/2. 

Outside  the  weakly  interacting  regime,  the  assumption  that  the  density  shift  scales  mostly  linearly  with  the  particle  number  breaks  down. 
One  important  consequence  of  this  observation  is  that  the  behavior  of  the  N  =  Nt  particle  system  cannot  be  reproduced  by  solving  the  N  =  2 
system  and  then  rescaling  the  interactions,  as  claimed  in  Ref.  [17],  i.e.  replacing  xa  and  in  the  N  —  2  solutions  by  Xif  — ►  Xfl (TVr  —  1) 
and  — y  J^Nt/2.  The  relevant  role  of  non-linear  effects  is  illustrated  in  Fig.  9.  In  the  top  panel  we  compare  the  behavior  of  three 

systems  with  different  N  =  2, 6, 12.  We  fix  the  interactions  such  that  Xa{N  —  1)  are  all  equal  as  well  as  all  J^N .  While  the  curves  N  =  6 
and  12  are  very  close  to  each  other,  they  do  not  agree  with  the  N  =  2  solution. 

Note  also  that  for  the  many-body  system,  N  ($>  2,  the  part  of  the  shift  sensitive  to  the  second  pulse  area  is  mainly  proportional  to 
(up  to  1/N  corrections)  and  not  to  2J^  +  Xfl  sin2(#i),  as  it  is  the  case  for  N  =  2.  This  phenomenon  is  visible  in  the  bottom  panel  of  Fig.  9, 
which  shows  the  shift  vs  excitation  fraction  for  three  different  second  pulse  areas:  (0.25,  0.5,  and  0.75)7r.  At  Jx  =  0  and  0"  =  n/2  (close 
to  zero  crossing  for  the  parameters  used  in  the  plot)  the  density  shift  only  exhibits  a  significant  dependence  on  the  second  pulse  area  when 
N  =  2. 
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FIG.  9.  (Color  Online)  Shift  vs  excitation  fraction  for  different  atom  numbers:  N  =  2  (Black  open  squares),  6  (no  symbol),  and  12  (red 
circles).  The  convention  for  the  second  pulse  area  is:  :  0.25-7T  (solid),  0.57T  (dashed  line)  and  0.757T  (short-dashed  line).  We  set  Cr  =  0. 

The  values  of  the  other  interactions  were  chosen  such  that  xj(]V  —  1)  are  all  equal  as  well  as  JrN.  For  N  =  2,  we  set  \r  =  1  and  J r  =  0 
and  —0.6,  for  the  top  and  bottom  panels  respectively.  Note  that  for  N  >  2  and  Jr  =  0,  the  density  shift  is  quite  insensitive  to  the  second 
pulse  area  as  predicted  by  theory. 


B.  TWA  and  Mean-field  solutions 


We  will  now  use  the  TWA  to  go  beyond  the  collective  approximation  that  was  required  to  derive  the  equations  of  motion  in  prior  sections. 
The  TWA  allows  us  to  include  different  types  of  mechanisms,  which  bring  the  system  out  of  the  collective  Dicke  manifold,  such  as  single¬ 
particle  excitation  inhomogeneities  and  non-collective  elastic  and  inelastic  two-body  collisions.  We  can  also  easily  include  one-body  losses 
due  to,  for  example,  background  gas  collisions.  The  latter  are  characterized  by  the  jump  operators  ^/T jacnja,  which  annihilate  the  atom  j  in 
spin  state  a. 

For  the  mean-field  treatment,  we  assume  the  reduced  density  matrix  of  the  many-body  system  factorizes  as  described  in  Eq.  (47).  With  the 
aim  of  dealing  with  long-range  interactions,  it  is  better  to  write  the  equations  of  motion  of  the  Fourier  transformed  quantities: 


N  2  ,fc 

paa{k)  =  y  pj^a 

(58) 

3=1 

N 

Peg(k)  =  Pi,g  =  P*ge(k) 

(59) 

3=1 

The  equations  of  motion  derived  in  this  way  are  general,  and  in  section  VII,  will  be  applied  to  model  the  Yb  clock  dynamics.  These 
equations  reduce  to  Eqs.  (49-51)  when  the  interaction  parameters  are  collective,  i.e.  have  only  k  =  k!  =  0  components.  In  their  most  general 
form,  the  equations  of  motion  are  nontrivial,  and  we  present  them  in  Appendix  4.  Here  instead  we  proceed  to  analyze  simpler  relevant  cases, 
which  will  allow  us  to  test  the  validity  of  their  corresponding  TWA  dynamics  (obtained  after  sampling  the  solutions  of  the  mean-field  equations 
with  the  Wigner  distribution).  We  proceed  by  considering  two  cases,  one  in  which  the  inhomogeneity  is  caused  by  interactions  and  the  other  in 
which  it  is  caused  by  single-particle  terms.  Since  analytic  solutions  are  only  available  when  one  restricts  the  dynamics  to  a  pure  Hamiltonian 
evolution,  in  the  reminder  of  this  section  we  will  neglect  single-particle  and  two-body  losses. 
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FIG.  10.  (Color  Online)  Ramsey-fringe-contrast  decay  for  the  case  of  power-law  Ising  interactions  in  a  2D  square  array  composed  of  121 
atoms  (See  text).  The  first  pulse  is  set  to  tr/2.  The  solid  lines  correspond  to  the  exact  solution  obtained  from  Ref.  [52],  The  dashed  lines  are 
the  TWA  solution.  From  top  to  bottom  the  curves  correspond  to  a  =  1,  0.5,  and  0.2. 
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FIG.  11.  (Color  online)  Top  panels:  Comparisons  of  the  density  shift  vs  excitation  fraction  predicted  by  the  collective  TWA  and  the  generic 
TWA  at  t  =  50  ms.  The  generic  TWA  accounts  for  the  mode  dependence  of  p-wave  interactions.  Here  zero  excitation  inhomogeneity  is 
assumed  (Afl  =  0).  A  thermal  average  over  initially  populated  modes  is  performed.  The  left  panel  is  for  a  fixed  number  of  atoms,  N  =  20.  In 
the  right  panel,  we  have  done  an  average  over  the  atom  number  in  the  pancakes  using  a  Poissonian  distribution  with  N  ~  9  (see  Sec.  VII  A). 
Here  vee  =  v" ,  thus  Cn^.n  ,  =  0.  Bottom  panels:  Comparison  of  the  contrast  decay  vs  time  predicted  by  the  collective  TWA  and  generic 
TWA  for  different  first  pulse  areas  (indicated  in  the  plot).  Again,  the  left  panels  is  for  a  fixed  number  of  atoms,  N  =  20,  and  the  right  panel  is 
computed  after  performing  an  atom  number  average. 


1.  Non-collective  interactions 


•  Ising  Case 

The  Ising  case  corresponds  to  .  n  ,  =  0.  An  exact  solution  for  the  Ising  dynamics  exists  for  an  initial  product  state  and  for  arbitrary 
spin  coupling  constants  [31,  52,  60-62].  This  analytic  solution  allows  us  to  benchmark  the  validity  of  the  TWA  beyond  the  collective  regime, 
assuming  zero  transverse  interactions  Jr| .  n  ,  =  0  and  no  losses.  We  will  also  assume,  for  simplicity,  no  excitation  inhomogeneity  during 
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FIG.  12.  (Color  online)  Left  panel:  Comparison  of  the  density  shift  vs  excitation  fraction  predicted  by  the  analytic  perturbative  treatment  (solid 
blue  line),  the  TWA  (filled  blue  dots)  and  the  mean-field  (red  squares)  with  excitation  inhomogeneity.  We  assumed  collective  interactions,  set 
Afi/f2  =  0.25  and  fixed  N  =  30.  The  density  shift  is  computed  at  =  0.04.  Here  vee  =  v",  thus  Cr  =  0.  Right  panel:  Ramsey 
contrast  decay  vs  time  including  excitation  inhomogeneity  predicted  by  the  analytic  perturbative  treatment  (solid  blue  line),  the  TWA  (filled 
blue  dots),  and  the  mean-field  (red  squares).  The  pulse  area  is  set  to  n/2. 


pulses.  In  this  case,  at  the  mean-field  level,  one  obtains: 
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Here  /  is  the  identity  matrix.  (Sq ,  Sq,  Sq}  are  the  components  of  the  collective  Bloch  vector  after  the  first  pulse.  In  Fig.  10,  we  show  the 
comparisons  for  a  system  of  N  =  121  =  11  x  11  atoms  assuming  that  modes  n  form  a  2D  square  lattice  array.  We  set  Xftj  ,fl  .,  1  / 1  n j  —  ny '  |  “ 

and  vary  a  =  0.2,  0.5, 1.  Note  this  choice  of  Xn,  is  done  just  for  convenience  without  any  direct  link  with  the  mode  dependence  of  Xr;j  ,h 
when  it  describes  actual  s  or  p-  wave  interactions.  The  physically  relevant  case  for  the  clock  will  be  described  below.  For  a  —  0.2, 0.5,  the 
exact  and  TWA  solutions  agree  well.  However,  for  a  =  1,  deviations  start  to  become  important.  An  approximate  criterion  for  the  validity  of 
the  TWA  for  a  translationally  invariant  system  is  that,  if  one  defines  [31]  Em  =  JW  xSJ,  Hj  >  then  E?  «  N E2.  Note  that,  in  this  analysis,  we 
have  not  added  any  quantum  vacuum  noise  to  the  initially  empty  k  ^  0  modes.  However,  we  have  confirmed  that  its  addition  only  slightly 
changes  the  conclusions. 

•  Pure  p- wave  interactions 

We  now  take  advantage  of  the  TWA  to  benchmark  the  validity  of  the  collective  model  when  the  interaction  inhomogeneity  comes  from  the 
mode  dependence  of  the  p-wave  interactions.  In  this  case,  for  simplicity  we  also  assume  no  excitation  inhomogeneity  during  the  pulses  and 
assume  no  single  particle  or  two-body  losses.  We  also  set  s-wave  interactions  to  zero  and  perform  a  thermal  average  over  populated  modes 
assuming  a  Boltzmann  distribution  (see  Sec.  VII  A).  In  Fig.  11,  we  show  the  density  shift  and  the  Ramsey  fringe  contrast  decay  predicted 
by  the  TWA  solution,  both  with  and  without  interaction  inhomogeneity.  While  the  density  shift  agrees  fairly  well  between  the  collective  and 
non-collective  solutions,  the  contrast  shown  for  the  high  (9 1  =  0.87r)  initial  pulse  area  has  a  slower  decay  for  the  case  of  collective  interactions. 
On  the  contrary,  at  initial  pulse  areas  close  to  n/2,  the  agreement  between  the  two  solutions  is  fair,  and  the  collective  model  exhibits  just  a 
slightly  faster  decay  of  the  contrast.  The  disagreement  of  the  contrast  at  high  excitation  fraction  is  because,  at  these  conditions  (a  similar  effect 
happens  at  small  tipping  angles),  the  contrast  reduction  comes  primarily  from  simple  dephasing  due  to  the  different  precession  rates  of  the 
Bloch  vector  for  different  realizations  of  the  interaction  parameters.  Recall  that  the  dynamics  induced  by  mean-field  interactions  correspond  to 
a  precession  around  an  effective  magnetic  field  whose  strength  depends  on  atom  number,  interaction  strength,  and  excitation  fraction,  Eq.  (39). 
For  the  collective  model  that  type  of  dephasing  is  much  weaker.  On  the  contrary,  if  Cr.  ^ =  0,  as  it  is  in  this  example,  the  effective  magnetic 
field  vanishes  at  9 1  =  n/2.  Thus,  for  9 1  close  to  n/2,  dephasing  is  minimal  and  contrast  decay  is  dominated  by  many-body  effects,  which  are 
fairly  well-captured  by  the  collective  model. 

Despite  the  disagreement  of  the  contrast  at  low/high  excitation  fractions  and  at  zero  s-wave  interactions,  the  collective  model  remains  very 
useful  for  describing  experimental  observations  during  the  interrogation  of  the  clock.  The  reason  is  that,  instead  of  a  single  pancake  with  fixed 
atom  number,  in  the  experiment,  there  is  an  array  of  pancakes  populated  with  a  Poissonian  distribution  (as  explained  below).  If  one  takes  into 
account  the  Poissonian  average,  one  recovers  a  good  agreement  between  the  collective  and  non-collective  models  even  at  low/high  excitation 
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fractions,  as  shown  in  Fig.  1 1.  In  this  case,  dephasing  due  to  the  different  atom  number  in  each  lattice  pancake  dominates  the  contrast  decay  at 
low/high  excitation  fractions.  When  the  Poissonian  average  is  included,  the  cases  with  low/high  pulse  areas  are  the  ones  that  exhibit  the  fastest 
Ramsey  fringe  contrast  decay.  Note  that,  for  Fig.  11,  the  mean-atom  number  per  the  pancake  used  for  the  plots  with  a  Poissonian  distribution 
was  N  ~  9.  On  the  contrary,  for  the  fixed-atom-number  plots  we  set  N  =  20.  The  difference  in  atom  number  is  the  only  reason  why  the 
overall  timescale  for  the  contrast  decay  in  the  Poissonian-averaged  plots  is  slower  than  in  the  fixed-atom-number  ones. 

We  also  expect  that  finite  |  J~.  H  / 1  >  0  interactions  (induced  by  s-wave  collisions)  will  enhance  the  gap  between  the  collective  Dicke 
manifold,  S  =  N/2,  and  the  other  manifolds.  A  larger  gap  will  further  suppress  population  leakage  outside  the  former  (gap  protection)[48] 
and  will  improve  the  collective  approximation.  The  validity  of  the  gap  protection  idea  is  explicitly  shown  in  Appendix  5,  where  we  use  a 
box  potential  instead  of  a  harmonic  potential  to  make  the  mode  dependence  of  the  p-wave  interactions  even  stronger.  There  we  show  that  the 
dynamics  converge  to  the  collective  dynamics  as  the  mean  value  of  |  Jrj.  -  J  is  increased. 


2.  Excitation  induced  inhomogeneity 


In  the  presence  of  finite  ip,  Eq.  (6),  the  gap  protection  does  not  help  to  prevent  leakage  of  population  outside  the  Dicke  manifold.  This  is 
because  the  transfer  is  made  during  short  pulses  where  interactions  play  no  role  and  single-particle  excitation  inhomogeneity  is  the  leading 
transfer  mechanism.  In  Fig.  12,  we  compare  the  analytic  solution  found  in  Sec.  VIA  using  perturbation  theory  in  Afi/D  (solid  blue  line), 
the  TWA  (filled  blue  circles),  and  a  pure  mean-field  with  excitation  inhomogeneity  (red  squares).  Fig.  12  assumes  collective  interactions 
and  equal  first  and  second  pulse  areas.  For  comparison  purposes,  we  also  show  the  analytic  solution  for  the  density  shift  in  the  case  of  zero 
inhomogeneity  (dashed  purple  line).  Note  that  excitation  inhomogeneity  manifests  itself  in  the  density  shift  mainly  at  high  excitation  fractions 
where  it  gives  rise  to  some  curvature.  The  figure  also  shows  that  the  density  shift  is  well  captured  by  both  the  TWA  and  the  mean-field.  On 
the  contrary,  the  Ramsey  fringe  contrast  decay  is  only  captured  by  the  TWA  and  analytic  solutions,  which  agree  well  with  each  other.  The 
mean-field  does  not  predict  any  contrast  decay.  The  pulse  area  for  the  Ramsey  fringe  contrast  shown  in  Fig.  12  was  set  to  be  n/2  . 


VII.  COMPARISONS  WITH  OPTICAL  LATTICE  CLOCK  EXPERIMENTS 


After  developing  and  benchmarking  in  great  detail  the  theoretical  formalism,  we  now  proceed  to  apply  it  to  describe  experimental  mea¬ 
surements  performed  in  the  8'Sr  JILA  optical  lattice  clock  [11,  22]  and  in  the  NIST  171  Yb  [19,  20]  optical  lattice  clock.  One  of  the  main 
differences  between  the  two  clocks  is  their  temperature.  While  the  JILA  clock  operates  typically  at  temperatures  between  (1  —  4)/iK  and 
with  a  typical  excitation  inhomogeneity  AS2/S2  <  0.1,  the  171  Yb  clock  typically  operates  at  temperatures  at  or  somewhat  above  10/rK  with 
corresponding  excitation  inhomogeneity  Af 1/Cl  ~  0.25  —  0.3.  These  differences  have  significant  consequences  for  the  proper  theory  model 
required  to  describe  the  experiment.  While  the  pure  collective  model  is  sufficient  to  describe  the  many-body  dynamics  in  the  8'  Sr  clock,  as 
explicitly  shown  below,  the  inclusion  of  excitation  inhomogeneity  is  needed  to  reproduce  the  1 7 1  Yb  measurements.  We  first  start  by  describing 
in  more  detail  the  experimental  set-up.  which  is  common  for  the  two  clocks,  and  then  present  the  experimental  measurements  and  comparisons 
with  theory. 


A.  Spatial  averaging  and  finite  temperature 

In  the  experiments,  a  one  dimensional  lattice  is  used  to  confine  the  atoms  creating  an  array  of  2D  disk-shaped  arrays  or  pancakes  (see 
Fig.  2b).  The  atom  number  distribution  along  the  array  is  determined  by  assuming  the  atoms  are  loaded  in  the  MOT  (magneto-optical  trap) 
with  a  Gaussian  density  profile,  and  that  this  profile  is  kept  the  same  but  simply  partitioned  into  lattice  sites  when  the  ID  lattice  is  turned  on. 
In  this  case  p(l)  =  NtZ  is  taken  as  the  density  at  pancake  l  after  integrating  along  the  other  two  directions.  Here  <jm  is 

the  1/e  MOT  radius,  which  is  ~  30/mi  for  8'  Sr  clock  and  ~  130/rm  for  the  171  Yb  clock,  a  ~  0.4  pm  is  the  lattice  spacing  (approximately 
the  same  for  both  experiments)  and  Nt  the  total  atom  number.  The  atom  number  fluctuates  from  realization  to  realization  according  to  a 
Poissonian  distribution,  and  the  probability  of  having  N  >  0  atoms  at  lattice  site  l  is  7Z(l,  N)  =  Ae~p^  with  A  a  normalization 

constant.  Any  experimental  observable  should  be  then  calculated  as  (DexP  =  J2N  0(N)1Z(N)  with  1Z(N)  =  71(1,  N).  In  Fig.  13,  we 
show  an  atom  number  distribution  for  a  typical  number  of  Nt  ~  26000  atoms  in  the  Yb  clock.  The  Sr  clock  operates  with  a  smaller  atom 
number,  Nt  =  5000,  but,  due  to  different  trapping  conditions,  the  average  distribution  of  atoms  in  both  experiments  is  similar,  as  can  be  seen 
in  Fig.  13. 

In  most  of  the  calculations  performed  above,  we  fixed  the  set  of  vibrational  modes  to  be  n.  However,  to  compare  with  experiments,  we  need 
to  perform  a  thermal  average.  For  /iK  temperatures,  at  which  the  clocks  operate,  we  can  assume  the  modes  are  initially  populated  according 
to  a  Maxwell-Boltzmann  distribution,  computed  using  single-particle  energies.  Although  so  far  we  have  assumed  that  atoms  are  frozen  along 
the  axial  direction,  this  assumption  is  not  fully  correct  for  the  Yb  lattice  clock  operated  at  ~  10/iK  with  a  typical  axial  trapping  frequency  of 
vz  =  75  —  80  KHz.  In  this  case,  the  mean  axial  mode  occupation  number  fiz  ~  2.5,  and  to  compare  with  the  Yb  experiment,  we  perform  a 
thermal  average  over  excited-optical-lattice  bands. 
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Atom  number  per  pancake 


FIG.  13.  (Color  Online)  Atom  number  distribution:  The  solid  blue  distribution  is  for  the  Yb  ID  optical  lattice  clock  at  a  typical  atom  number 
Nt  =  26000;  the  dashed  purple  distribution  is  for  the  Sr  ID  optical  lattice  clock  at  typical  atom  number  Nt  =  5000. 


Excitation  fraction 


FIG.  14.  (Color  Online)  Density  shift  in  Ramsey  spectroscopy  vs  excitation  fractions  measured  in  the  Sr  clock:  The  filled  symbols  show 
experimental  data.  The  dashed(blue)  line  is  the  prediction  from  the  collective  TWA  model  and  the  solid  (green)  line  the  collective  mean-field 
model.  The  experimental  data  was  normalized  to  a  typical  time-averaged  total  atom  number  of  1000  atoms.  The  excitation  fraction  was  also 
time-averaged.  See  Ref.  [22]  and  text  for  more  information. 


B.  Ramsey  Spectroscopy 

The  experiments  performed  Ramsey  spectroscopy  measurements  and  varied  the  atom  number,  pulse  area,  and  the  initial  state  of  the  atoms. 
The  results  are  summarized  in  Figs.  14-16  for  the  Sr  clock  and  Figs.  17-18  for  the  Yb  clock.  In  those  plots,  the  theoretical  predictions  are  also 
shown. 


1.  Sr  optical  lattice  clock 


Let  us  first  discuss  the  Sr  clock,  which  admits  a  simpler  theoretical  treatment  based  on  a  pure  collective  model.  Figure  14  shows  the  measured 
density  shift  vs  excitation  fraction.  The  excitation  fraction  was  changed  by  varying  the  first  pulse  area.  In  the  presence  of  two-body  losses, 
the  excitation  fraction  is  not  constant  during  the  dark  time,  and,  to  simplify  the  comparisons  with  theory,  a  time-averaged  excitation  fraction 
was  used  for  the  density  shift  plot.  The  latter  was  extracted  in  the  experiment  by  performing  a  set  of  independent  measurements  periodically 
inserted  into  the  clock  sequence  [22].  The  duration  of  the  dark  time,  r,  was  set  to  80  ms,  and  the  final  pulse  area  to  tt/2.  The  typical  axial 
and  radial  sample  temperature,  measured  by  sideband  spectroscopy  and  Doppler  spectroscopy  respectively,  were  between,  Tz  ~  1  —  2/rK  and 
Tii  ~  2  —  4/rK;  the  trapping  frequencies  were  set  to  vz  =  80  kHz  and  vr  ~  450  Hz;  and  the  typical  misalignment  angles  were  small  enough 
that  excitation  inhomogeneity  was  negligible.  Under  these  operating  conditions,  the  collective  model,  after  performing  an  average  over  the 
atom  number  distribution,  is  expected  to  work  fairly  well.  This  expectation  is  confirmed  by  the  linear  dependence  of  the  shift  at  high  excitation 
fraction  and  the  good  agrement  between  the  collective  TWA  predictions  and  the  experimental  data.  Since  only  the  g  interaction  parameters  are 
known,  bgg  ~  76  as  (as  the  Bohr  radius)[ll],  we  used  the  density  shift  data  to  determine  the  parameters  Ve’9  —  V 9,9  =  —  2-7T  x  0.27  Hz 
and  Ve'e  —  V 9,9  =  0AVe’3 .  Note  that  these  parameters  are  almost  temperature  independent. 

Figure  15  shows  the  atom  number  decay  and  the  normalized  Ramsey  fringe  contrast.  By  normalized  we  mean  that  the  contrast  was  rescaled 
by  the  corresponding  atom  number  at  a  given  time.  In  the  experiment,  the  Ramsey  fringe  contrast  was  extracted  in  a  manner  that  is  insensitive 
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FIG.  15.  (Color  Online)  Atom  number  decay  and  Ramsey  fringe  contrast  vs  dark  time  measured  in  the  Sr  clock:  The  filled  symbols  show 
experimental  data,  the  solid  lines  are  the  predictions  from  the  collective  TWA  model,  and  the  dashed  lines  are  the  predictions  from  the  collective 
mean-field  model.  See  text  for  more  information.  An  overall  scaling  factor  on  the  order  of  ~  0.9  —  1  in  the  contrast  plots,  that  accounted  for 
the  experimental  dead-time  between  the  end  of  the  second  pulse  and  the  population  measurements  was  used  in  the  theory  curves  to  match  the 
contrast  at  zero  dark  time.  Statistical  error  bars  for  experimental  data  are  comparable  to  symbol’s  size.  See  also  Ref.  [22], 


to  the  frequency  noise  of  the  local  oscillator.  For  details,  see  Ref.  [22,  23] .  Four  different  first  pulse  areas  were  used  (the  second  pulse  area  was 
always  n/2)  and  are  indicated  with  different  symbols  in  the  figure.  The  Ramsey  pulses  were  always  shorter  than  6  ms.  Atom  number  decay 
data  was  used  to  determine  Te,e  =  27r  x  0.1  Hz.  Fe’3  was  consistent  with  zero  in  agreement  with  measurements  reported  in  Ref.  [21].  There 
the  single  atom  decay  rates  were  also  measured  to  be  Fe  ~  T9  =  27r  x  0.01  Hz,  and  we  used  those  values.  These  parameters  together  with 
the  values  of  Ve'3  —  V3’3  and  Ve'e  —  V3’3  extracted  from  the  density  shift  were  used  to  compute  the  Ramsey  contrast  decay  curves. 

We  added  a  cubic  term  to  account  for  the  virtual  populations  of  the  off-resonant  modes  as  explained  in  Sec.V.  Since  the  coupling  constants 
af  2j3  cannot  be  trivially  computed,  we  used  instead  the  quantities  Xp’p ^  ~  which  give  rise  to  the  relevant  cubic  contributions  -  as  fitting 
parameters  (see  Appendix  2)  and  set  a^g  =  0.  The  values  of  Xpp]  utilized  generated  an  additional  term  Ar(Ne) 3  r(Ng)r  with 

Ao  =  27r  x  0.0016  Hz,  Ai  =  —2n  x  0.0058  Hz,  A2  =  2n  x  0.003  Hz,  A3  =  27r  x  0.0011  Hz.  Although  the  cubic  corrections  are  small 
and  played  no  role  in  the  density  shift,  they  improved  the  agreement  between  theory  and  experiment.  To  emphasize  this  point  in  Fig.  16,  we 
explicitly  compare  the  theory  predictions  of  the  contrast  with  and  without  cubic  terms. 

In  addition  to  the  collective  TWA  solutions,  we  also  show  the  collective  mean-field  predictions.  Both  were  computed  after  performing  an 
average  over  the  atom  number  distribution.  Excellent  agreement  is  observed  between  the  collective  TWA  and  the  experimental  data  [63].  The 
mean-field  model  reproduces  the  density  shift  measurements  and  is  able  to  capture  the  fast  contrast  decay  at  low  excitation  fractions  due  to 
the  dephasing  induced  by  the  different  precession  rates  of  the  Bloch  vector  across  different  pancakes.  Remember  that  the  effective  magnetic 
field  induced  by  interactions  is  directly  proportional  to  the  density  shift  and  depends  on  the  pulse  area  and  atom  number.  At  low  excitation 
fraction  the  effective  magnetic  field  (density  shift)  is  maximum  as  shown  in  Fig.  14.  On  the  contrary,  for  the  intermediate  excitation  fraction 
data  points,  the  mean-field  solution  fails  to  reproduce  the  contrast  decay  observed  experimentally.  At  those  excitation  fractions,  the  effective 
magnetic  field  (density  shift)  is  almost  zero  -see  Fig.  14  -  and  the  contrast  decay  is  dominated  by  many-body  correlations  not  included  in  the 
mean-field  model. 


2.  Yb  optical  lattice  clock 

Now  we  proceed  to  discuss  the  Yb  clock  experiment.  Figure  17  shows  the  measured  density  shift  vs  excitation  fraction  seen  in  the  Yb 
clock  [19,  20].  Two  types  of  initial  conditions  were  considered.  For  the  red  and  blue  points,  initially  all  atoms  were  prepared  in  the  g  state 
and  then  optically  excited  to  the  e  state.  The  blue  data  (solid  circles)  was  taken  using  a  vertical  optical  lattice  and  the  red  data  (solid  squares) 
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FIG.  16.  (Color  Online)  Ramsey  contrast  decay  obtained  from  the  the  collective  TWA  with  cubic  terms  (solid  lines)  and  without  them  (dashed 
lines).  The  symbols  are  experimental  data  taken  at  different  first  pulse  areas.  Statistical  error  bars  for  experimental  data  are  comparable  to 
symbol’s  size.  See  also  Ref.  [22], 
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FIG.  17.  (Color  Online)  Density  Shift  vs  excitation  fractions  measured  in  the  Yb  clock:  The  filled  symbols  show  experimental  data  and  the 
open  symbols  joined  by  dashed  lines  show  predictions  from  the  non-collective  TWA  model  (left)  and  the  non-collective  mean-  field  model 
(right).  The  blue  circles  and  red  squares  were  taken  using  a  vertical  and  horizontal  ID  lattice,  respectively,  and  used  a  g  — >  e  interrogation. 
The  yellow  diamonds  used  a  vertical  lattice  but  interrogated  the  e  — >  g  transition.  The  experimental  data  measured  the  shift  by  normalizing 
it  to  a  typical  time-averaged  total  atom  number  of  (Nt)t  =  25000  atoms.  The  excitation  fraction  was  also  time-averaged.  See  text  for  more 
information. 


using  an  horizontal  optical  lattice.  In  the  vertical  lattice,  tunneling  between  neighboring  lattice  sites  is  expected  to  be  substantially  suppressed 
by  gravity  [19,  64],  That  suppression  is  not  present  in  the  horizontal  setup.  Based  on  that  and  the  fact  that  there  is  no  substantial  difference 
between  the  two  sets  of  data,  we  conclude  that,  at  least  for  the  density  shift,  tunneling  can  be  neglected.  The  yellow  solid  diamonds  show  the 
density  shift  measured  for  the  initial  condition,  in  which  all  the  atoms  were  prepared  in  the  e  state  and  then  optically  deexcited  to  the  g  state. 
For  all  the  cases,  the  Ramsey  pulse  time  was  fa  =  fa  =  1  ms,  the  dark  time  was  r  =  80  ms,  the  temperature  was  T  =  lOg  K,  the  trapping 
frequencies  were  vz  =  75  kHz  and  vr  ~  500  Hz  and  Afi/fi  ~  0.25.  In  171  Yb,  both  e  —  e  and  e  —  g  losses,  as  well  as  single  particle  e 
and  g  losses,  play  a  non-negligible  role  in  the  quantum  dynamics  [20].  Consequently,  both  the  excitation  fraction  and  atom  number  vary  with 
time  during  the  spectroscopy.  The  data  presented  here  accounts  for  those  effects  by  using  time-averaged  excitation  fraction  and  a  density  shift 
normalized  to  a  time-averaged  atom  number.  The  different  behavior  observed  between  the  e  —  g  and  g  —  e  interrogation  measurements  is 
mainly  due  to  excitation  inhomogeneity,  which  depends  on  pulse  area. 

Figure  18  shows  the  atom  number  decay  and  normalized  Ramsey  fringe  contrast.  Three  different  pulse  areas  were  used  (both  the  first  and 
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FIG.  18.  (Color  Online)  Atom  number  and  Ramsey  fringe  contrast  vs  dark  time  measured  in  the  Yb  clock:  The  filled  symbols  show 
experimental  data,  the  solid  lines  are  the  predictions  from  the  non-collective  TWA  model,  and  the  dashed  lines  are  the  predictions  from  the 
non-collective  mean-field  model.  See  text  for  more  information.  To  account  for  the  additional  25  ms  dead-time  between  the  end  of  the  second 
pulse  and  the  population  measurements,  we  introduced  an  overall  scaling  factor  ~  0.8  —  0.9  to  match  the  experimentally  measured  contrast  at 
zero  dark  time.  Experimental  error  bars  are  comparable  to  symbol’s  size. 


second  pulse  areas  were  the  same)  as  indicated  with  different  symbols  in  the  figure.  For  each  pulse  area,  two  measurements  were  done,  one 
at  a  high  and  the  other  at  a  low  atom  number.  These  measurements  were  carried  out  at  similar  but  slightly  different  trapping  conditions  and 
temperature  as  compared  to  the  ones  used  for  the  density  shift:  vz  =  80  kHz,  i/r  ~  550  Hz,  T  >  10/rK,  and  t\  =  t2  =  5  ms. 

The  most  general  TWA,  which  accounts  for  excitation  inhomogeneity,  non-collective  interactions,  two-body  and  single-particle  losses,  as 
well  as  interaction  effects  and  virtual  population  of  motional  states  (cubic  terms  in  the  spin  Hamiltonian),  was  used  to  model  the  experiment. 
Both  a  thermal  average  and  an  average  over  the  atom  number  distribution  in  the  pancakes  were  performed.  We  also  compared  the  predictions 
of  the  pure  mean-field  model,  which  accounts  for  similar  types  of  effects  but  without  performing  the  sampling  over  the  Wigner  distribution. 
Only  the  g  interaction  parameters  are  known  in  Yb  as  well,  bgg  «  0,  and  thus  we  used  ve’e,  ve,B ,  Ye’e,  ueg,  X'pp'1  (see  Appendix  2)  as 
fitting  parameters.  Note  that,  since  we  are  not  using  the  collective  model,  we  quote  the  v“,/3  values  and  not  the  values.  The  cubic 
Hamiltonian  was  the  only  term  projected  onto  the  collective  manifold  and  included  as  a  term  ar[N^~r {k  =  0 )Ng(k  =  0)].  Here 

P  =  1/iV2  Pn.  ,n .  ar  are  functions  of  Xpp>Xpp-  From  Ref.  [20],  we  set  7®’®  =  37e,e/5  and  Te  ~  T®  =  27r  x  0.3  Hz.  The 
theoretical  curves  shown  in  Figs.  17  and  18  used  the  following  set  of  parameters:  ve’9  =  —2n  x  3.8  Hz,  ve,e  =  — 0.2ve’®,7e,e  =  2-k  x  2.8 
Hz,  ueg  =  0,  ao  =  03  =  2n  x  0.3  Hz  and  ai  =  (13/2  =  —  2ao/3.  For  the  Ramsey  contrast  measurements,  we  found  that  the  5  ms  pulses 
used  in  the  Yb  experiment  were  not  short  enough  to  make  interaction  effects  negligible  during  the  pulses,  and,  because  of  that,  we  included 
interactions  during  those  in  the  numerical  simulations.  We  added  in  addition  a  single-particle  dephasing  term  acting  on  peg  in  all  cases,  e  UpT 
with  7 dP  =  2  Hz  to  capture  the  contrast  decay  observed  for  the  low  density  data  since,  at  those  densities,  interactions  effects  could  not  be  the 
cause  of  contrast  decay. 

While  there  is  a  substantial  number  of  fitting  parameters,  there  is  also  a  variety  of  experimental  data  curves.  Note  also  that  •ye,e  is  uniquely 
determined  by  the  atom  loss  curves  and  the  ratio  ve,e/ve’9  by  the  density  shift  measurements.  In  the  Sr  experiment,  the  collective  model  based 
on  the  TWA  was  able  to  fairly  reproduce  the  data  taken  under  different  conditions  using  the  same  set  of  two-body  interaction  parameters  and 
decay  rates.  For  the  Yb  experiment,  the  theory  could  reproduce  the  experimental  measurements  of  the  contrast  and  density  shift  using  one 
set  of  interactions  and  loss  parameters,  up  to  a  rescaling  by  a  factor  of  2.5  in  total  atom  number  when  computing  the  density  shift  curves. 
The  factor  of  2.5  difference  could  in  principle  come  from  the  different  experimental  conditions  used  for  the  contrast  decay  and  density  shift 
measurements,  although  most  likely  it  has  its  roots  in  additional  decoherence  mechanisms  not  accounted  for  by  the  model,  such  as  tunneling, 
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higher  vibrational  bands,  and  motional  decoherence.  Those  are  likely  to  affect  the  Ramsey  fringe  contrast  decay  more  than  the  density  shift 
and  could  be  responsible  for  the  overestimation  of  the  interaction  parameters  when  we  tried  to  fit  the  contrast  decay.  We  did  find  that,  if  we 
set  all  dr  —  0  (no  cubic  terms),  we  needed  a  re-scaling  factor  much  larger  than  2.5  to  fit  all  the  experimental  data  with  the  same  interaction 
parameters. 

While  both  the  TWA  and  the  mean-field  model  reproduced  fairly  well  the  density  shift,  the  mean-field  model  did  a  poorer  job  capturing 
the  contrast  decay  than  the  TWA.  For  Yb,  the  effective  magnetic  field  vanishes  at  pulse  areas  close  to  7r/2,  and  notably,  it  is  at  this  pulse  area 
when  the  deviation  between  the  mean-field  and  the  TWA  is  more  significant.  Similar  behavior  was  seen  when  modeling  the  Sr  clock  and  is 
consistent  with  the  expectation  that  when  the  effective  magnetic  field  induced  by  interactions  approaches  zero,  quantum  correlations  are  the 
dominant  cause  of  Ramsey  fringe  contrast  decay.  It  suggests  that  even  at  T  =  10/iK,  quantum  correlations  can  be  playing  a  non-negligible 
role  in  the  dynamics  of  the  interrogated  atoms. 


VIII.  CONCLUSIONS 


Here  we  have  derived  a  powerful  theoretical  formalism  capable  of  dealing  with  the  non-equilibrium  many-body  dynamics  of  open  spin 
models  with  long-range  interactions.  The  formalism  is  based  on  the  TWA,  which  we  have  applied  to  open  quantum  systems.  We  have 
benchmarked  the  accuracy  of  the  TWA  by  taking  advantage  of  existing  exact  solutions  as  well  as  by  comparisons  with  analytic  perturbative 
treatments  and  by  numerically  solving  the  master  equation.  We  have  applied  the  developed  formalism  to  describe  the  many-body  dynamics  of 
optical  lattice  clocks  during  Ramsey  spectroscopy  and  demonstrated  that  a  full  treatment  of  quantum  correlations  reproduces  best  the  observed 
dynamics,  better  than  a  pure  mean-field  treatment. 

The  formalism  developed  here  should  be  useful  for  the  description  of  a  broad  range  of  important  modern  quantum  systems,  including  for 
example,  trapped  ions,  neutral  atoms  in  optical  cavities,  quantum  dots,  and  nitrogen  vacancy  centers. 

Going  beyond  the  spin  model  formulation,  which  relies  on  the  assumption  of  frozen  motional  single  particle  states,  an  alternative  approach 
is  needed  to  deal  with  colder  samples  where  interactions  can  dominate  over  the  single-particle  energy.  In  this  parameter  regime  the  spin  model 
assumption  is  not  necessarily  valid.  One  possible  alternative  approach  includes  a  kinetic-theory  treatment  based  on  solving  a  Boltzmann 
equation.  That  type  of  treatment  has  been  shown  to  be  useful  for  describing  [65]  the  spin  segregation  observed  at  IILA  during  Ramsey 
spectroscopy  with  a  two-component  ultra-cold  bosonic  gas  [66]  and  for  describing  [67]  similar  segregation  effects  observed  at  Duke  in  a 
two-component  cold  trapped  Fermi  gas  [68]. 
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APPENDIX  1:  COLLECTIVE  MASTER  EQUATION 


In  this  Appendix,  we  present  the  details  behind  the  master-equation  treatment  of  Sec  IV  A.  Let  px  be  the  density  matrix  for  a  single  sector  of 
A f  particles.  As  explained  in  the  main  text,  since  there  is  no  coherence  between  the  sectors  with  different  particle-number,  the  master  equation 
can  be  solved  in  a  “block-diagonal  way”.  Specifically,  if  the  system  starts  with  N  particles,  we  need  to  solve  a  series  of  differential  equations 
for  each  of  the  subspaces  with  cascading  atom  numbers:  first  for  N  particles,  next  for  N  —  2  particles,  then  for  N  —  4  particles,  etc.  There  are 
(n)  =  n!(N-n)!  different  sectors  with  N  —  n  particles,  n  =  0,  2  . . .  N.  For  mode  independent  elastic  and  inelastic  interaction  parameters, 
in  each  of  the  sectors  with  Af  =  N  —  n  particles,  the  dynamics  are  restricted  to  the  collective  Dicke  states  |S  =  Af/2,  Mj v)  =  \Mu). 
Moreover,  each  of  the  sectors  behaves  identically.  If  we  assume  that  particles  in  pM  are  numbered  from  1  to  A/”  in  such  a  way  that  atoms 
Af  +  2  and  Af  +  1  are  the  ones  that  decay  as  one  goes  from  pM+2  to  pj^,  the  resulting  equations  are 
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for  0  <  Af  <  iV  —  2.  In  terms  of  spin  operators, 
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After  the  dynamics  of  the  various  density  matrix  elements  {M'm\pm\Mm)  are  known,  one  can  compute  any  observable.  For  example,  the 
number  of  e  atoms  is  given  by 
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where  Af  runs  over  Af  =  N,N  —  2,  TV  —  4, .... 


APPENDIX  2:  EFFECTIVE  HAMILTONIAN 


In  this  Appendix,  we  explain  in  detail  the  steps  that  lead  to  Eq.  (55)  in  Sec.  V.  As  explained  there,  the  states  I'I'J.)  can  be  divided  into  six 

k 

different  categories,  according  to  the  intermediate  state  participating  in  2nd-order  perturbation  theory.  Each  of  these  six  cases  gives  rise  to  a 
different  contribution  in  the  Hamiltonian;  we  will  now  present  these  contributions.  Let  define  ?  the  set  of  initially  unoccupied  modes. 


1.  Case  l^) 

k 

We  start  with  process  I,  where  in  the  intermediate  state,  one  atom  is  populating  a  mode  that  does  not  belong  to  the  initially  populated  manifold. 
We  obtain  terms  with  one,  two,  and  three  spin  operators,  i.e. 
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2.  Case 

k 

Next,  we  investigate  the  case  where  two  atoms  occupy  two  different  states  outside  the  initially  populated  modes.  We  arrive  at  terms  with  one 
and  two  spin  operators: 
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3.  Case 

k 

Now  we  consider  processes  of  type  2hlm  where  two  atoms  occupy  the  same  mode  outside  the  initial  configuration  space: 
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4.  Case  1 4^ ) 

k 

In  the  case  of  processes  of  type  Id,  two  atoms  occupy  the  same  mode  within  the  initial  configuration  space.  We  obtain  contributions  with  one, 
two,  and  three  spin  operators: 
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5.  Case  1^) 

k 

Next,  we  consider  the  case  when  two  modes  inside  the  initial  configuration  space  are  doubly  occupied.  We  obtain  contributions  with  one,  two 
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and  three  spin  operators: 
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6.  Case 
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Finally,  we  investigate  the  case  of  double  occupancy  of  a  mode  G  <?  and  one  mode  ^ 


rjlhld 

Ea 


E 


e« 


-En,  —  Pnh  +  -Eii,  —  -E'no 


-if' 


(lhld) 


(83) 


=  ^n3n,nhn,((Ve)2  - 

-  uegSl.nm{Sl.Sli  +  SI.SI,  +  Szn.S^) 

+  (vegfPnjn,nhnq(S^.S^  +  S9.Sl)  +  ((vee)2  -  (V9)2  +  K9)2)Pn3.n,nhn9^.^ 

+  Pn.n.n^n,,  [-«'>“  +  +  59 j <  +  +  S^S9 J  -  ((nee)2  +  (n")2)(^^  +  j] 

+  «e9Snj„,nfen,PnJnlnlln,  [(t>“  +  »W)(^^  +  -  S^h  -  S^S^J  -  2  V*9  (S^S^  -  S*  .  S*  j] 

_  9C„ee  _  ,,99i,,e9  P2  fO31  <?X  4_  0«  qz  i  O31  Cz  I  CK  OK  \ 

+  2(n6e  -  V39)UegPnjninhnqSnjninhnq(SnjSnhSni  +  -  S3J9J^.) 

-  2((vee)2  ~  (v")2)P2jnmS^S^Xh- 


(84) 


The  formulas  for  the  contributions  of  the  effective  Hamiltonian  can  be  simplified  significantly  if  we  consider  the  case  of  collective  interac¬ 
tions,  i.e. we  use  the  fact  that  the  coupling  constants  Pnj ninhnq  and  Snjninh nq  do  not  depend  significantly  on  i,  j ,  h,  and  q. 


Under  the  assumption  of  equal  coupling  constants,  it  follows  immediately  that  the  contributions  to  the  effective  Hamiltonian  stemming  from 
processes  of  type  Id  and  2d,  Eq.  (76)  and  Eq.  (81),  vanish  due  to  symmetry  considerations.  We  then  arrive  at  the  following  contributions: 
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Finally,  we  make  use  of  the  fact  that  within  the  collective  Dicke  manifold  the  relation  (S*)2  +  ( Sy )2  +  (Sz)2  =  y  (y  +  1)  holds. 
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l  h  1  (3) 

“2  =  2 Xpp 
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a2  =  Xpp 


((vey  -  (v33)2)2(N  —  1) 


es\2 


,,99\2', 


2({veey  -  2(vesy  +  (v"Y) 


i2hlm  r>/  \2  (5) 

b2  =  -2(ueg)  XsS’ 


a\hld  =  *(vee  -  v")  [x^piN2  +  3N-  2)(uee  +  v")  +  4JV(JV  -  4 )(x%ueg  -  x^V3)]  , 
a\hld  =  -(N  -  2)  [xpp((vee)2  +  2 (ve9Y  +  (v"Y  ~  2ve3(vee  +  v "))  +  2X<p)sueg(vee  +  2veg  +  u99)]  , 
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(100) 

(101) 
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(103) 

(104) 

(105) 

(106) 
(107) 


In  conclusion,  the  virtual  processes  give  rise  to  an  effective  Hamiltonian,  which  after  projection  onto  the  collective  Dicke  Manifold  and  up 
to  constants  of  motion  ( i.e .  S  ■  S  terms)  is  given  by: 


=aJSz+a^(SzY  +  aUSz)3, 

T  lh  ,  2h2m  ,  lhld 

a1  =  ax  +  , 

T  1  h  ,  2h2m  ,  lhld 

a2  =  a2  +  a2  +  ax  , 

T  lh  ,  lhld, 

a3  =  a2  +  or 


(108) 

(109) 

(HO) 

(111) 


APPENDIX  3:  ANALYTIC  SOLUTION  FOR  THE  CASE  OF  WEAK  EXCITATION  INHOMOGENEITY 


In  the  presence  of  excitation  inhomogeneity,  an  analytic  treatment  based  on  perturbation  theory  can  be  performed  when  dealing  with  purely 
unitary  evolution  (i.e.  neglecting  two-body  losses)  and  assuming  collective  two-body  interactions. 

To  accomplish  that,  one  writes  =  On  +  <5 fig,. ,  with  fig  =  y^.;  fig^  /N  the  mean  Rabi  frequency  and  treats  /Cl  as  a  perturbation 
parameter.  Note  that  in  this  case  the  Hamiltonian  evolution  during  the  pulse  can  be  written  as 


(u>:j  =eitn§ym, 

it(t)  =  (i  +  itd-  (tdy/2...), 


3=1 


%2-Kkj/N  Qy 


3  =  1 


(112) 

(113) 

(114) 
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Here  k  =  1  ,...N  -  1,  Sx’y'z 


,  Sflk 


1  v-'AT 
N  2sj=l 


■  2-i rj  fc  /'i 

e”  4  0(V,  and  Si 


=  5g.  ±*5H  .  Since 

rl.7  nj 


X)  ei2,rfej'/JV5J  |JV/2,  M)  =  A%m\N/2  -  1,  M  ±l,k), 

3  = 1 

„±  _  ,  AiV/2T.M)(jV/2TM-l)V/2 

^jv,m-±^  (AT-1)  J  ’ 


(115) 


then  0|7V/2,  M)  =  JZfcLi  <Sfifc(-A]v  M|iV/2  —  1,  M  +  1,  fc)  —  M|iV/2  —  1,  M  —  1,  k)).  This  equality  means  that  the  role  of  in¬ 
homogeneity  is  to  populate  the  S  =  N/2  —  1  states,  the  so  called  spin-wave  states.  These  states  during  the  free  evolution  acquire  a  phase 
—6 Mt  +  XflA/2T  ~  N  ji;r,  which  will  give  rise  to  an  additional  interaction  energy  shift  proportional  to  J not  present  in  the  fully  collective 
case.  5  =  6  —  Ca{N  —  1). 


Under  this  approximation,  one  obtains  that,  after  the  second  pulse,  the  number  of  excited  atoms,  2Ve,H  (ti ,  fe)  in  Ramsey  spectroscopy  is: 


lVe,  a(ti,t2)  =  y  +  (Sz)  =  N/2  +  cOs(6>ii)(^(0-)|5Rt(fi)eiT"S/fisRt(f2)5z5R(f2)e-<T"S/ftsR(f1)|V;(0-)) 

+  sin(0J)('0(O_)|3it(ti)e*T^  '/fi5Rt(f2)5a:K(f2)e_l'r^  ^fiSR(fi)|i/)(0_)}.  (116) 

Keeping  only  the  leading  order  correction,  which  can  be  shown  to  be  proportional  to  A .  /N  —  fl|,  i.e.  the  first-order  corrections 
vanish,  one  can  show  that  the  Ramsey  signal  is  given  by 


NeA(tiM)  =  —  A0iA02IUcos  (0i )  cos  (02)  ZN~2  cos[r(0  +  (N  -  2)C  +  0)] 

+  y  sin  (0i )  sin  (02)ZJV_1  cos  |r(<5  +  (N  -  l)C)j 

(N  -  1)2  cos  [t(<5  +  (N  —  3)C  +  0)]  +{N-  3)  cos  [r(<5  +  {N  —  3)C)] 


N 

Y 


Adi 


sin(0i)  sin (02 ) 2^ 


+  K, 


(117) 


The  parameters  j),  </>,  2,  W,  and  A'  are  given  by: 


tan(t)r)  =  tan  [^(ArJ±  +  x)tJ  sec(0i),  (118) 

tan (4>t)  =  tan  (2xt)  cos(0i),  (119) 

22  =  1  —  sin2(0i)  sin2  (2 \t)  ,  (120) 

W2  =  1  —  tan2(0i)  sin2  |^(ArJ±  +  x)rj  ,  (121) 

K  =  — - —  cos(0i)  cos (02)  (  1 - - 1-  j - —  sin(0i)  sin(02)A02A0i  cos (NJ±t).  (122) 


and  2,  </  given  by  Eqs.  (21)  and  (22).  For  simplicity,  we  have  removed  the  subscript  n  but  it  is  understood  that  all  the  interactions  and  Rabi 
frequencies  are  for  the  selected  modes  n  under  consideration.  For  N  =  1,  AQ,  =  0. 


APPENDIX  4:  GENERIC  MEAN-FIELD  EQUATIONS  OF  MOTION 


Here  we  derive  the  most  general  equations  of  motion  that  include  excitation  inhomogeneity,  non-collective  elastic  and  inelastic  two-body 
interactions,  and  single-particle  losses.  The  equations  are  obtained  by  assuming  that  the  density  matrix  can  be  factorized  according  to  Eq.  (47). 
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Here  we  have  introduced  the  following  quantities: 
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We  have  also  assumed  translationally  invariant  single  particle  decay  and  detuning  since  this  is  the  most  relevant  experimental  case.  However, 
those  can  be  made  to  be  inhomogeneous  straightforwardly.  At9  a  t  accounts  for  s-wave  e-g  inelastic  collisions.  We  have  not  introduced  s— 
losses  before  since  those  are  only  possible  when  the  dynamics  are  not  restricted  to  take  place  in  the  fully  symmetric  Dicke  manifold  [54], 

Note  that  for  a  coherent  state  ,  pee(k  =  0)pgg(k  —  0)  =  peg{k  =  0 )pge{k  =  0).  This  is  an  important  equality  to  keep  in  mind  to  reproduce 
the  terms  proportional  to  Te9  in  Eqs.  (49-51)  and  to  reproduce  the  absence  of  s-wave  losses  in  the  collective  limit. 
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APPENDIX  5:  GAP  PROTECTION 


o 


FIG.  19.  Effect  of  s-wave  interactions  on  the  validity  of  the  collective  approximation.  Comparison  of  the  exact  solution  (solid)  with  the 
collective  approximation  (dashed)  in  the  absence  of  s-wave  interactions  (top)  and  for  u/veg  =  15  (bottom). 


The  derivation  of  the  collective  spin  model  is  based  on  the  assumption  that  all  particles  interact  collectively  with  each  other  and  pulses  are 
homogeneous,  i.e.  the  coupling  constants  Pnjninhnq  and  Snjninhnq  are  taken  to  be  independent  of  j,  i,  h  and  q.  This  assumption  is  justified 
if  there  exists  a  large  energy  gap  that  prevents  transitions  between  the  S  =  N/2  and  S  =  N/2  —  1  sectors.  Such  a  gap  can  be  generated  if  the 
term  J^S  ■  S  is  large  compared  to  the  other  terms  in  the  Hamiltonian.  Recall  that  =  (VS9  —  US9)/ 2,  where  VS9  and  US9  are  functions 
of  the  scattering  lengths  for  p— and  s-wave  interactions,  respectively.  We  expect  that  the  assumption  of  a  collective  interaction  approximates 
well  the  full  dynamics  of  the  system  in  the  presence  of  a  large  energy  gap.  As  an  illustration,  we  present  in  Fig.  19  a  comparison  of  the  exact 
solution  with  the  collective  approximation  in  the  case  of  no  s-wave  interaction,  u  =  0,  and  a  ratio  of  u/vea  =  15.  In  this  case,  we  use  a 
boxed  potential  to  compute  the  mode-dependence  of  the  p-wave  interactions.  The  s-wave  interactions  in  the  boxed  potential  become  mode 
independent. 
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